Repositorio Institucional de Documentos

Universidad de Zaragoza Escuela de I Departamento de Ingeniería Electrónica y Comunicaciones Proyecto Fin de Carrera Ingeniería de Telecomunicación Especialidad en Comunicaciones SISTEMAS DE OPTIMIZACIÓN DIRECTA DE PARÁMETRO PROCESOS DE RADIOTERAPIA PARA EL TRATAMIENTO DEL CÁNCER Autor Director Zaragoza, Mayo de 2012 Universidad de Zaragoza Escuela de I ngeniería y Arquitectura Departamento de Ingeniería Electrónica y Comunicaciones Proyecto Fin de Carrera Ingeniería de Telecomunicación Especialidad en Comunicaciones SISTEMAS DE OPTIMIZACIÓN DIRECTA DE PARÁMETRO S DE MÁQUINA EN PROCESOS DE RADIOTERAPIA PARA EL TRATAMIENTO DEL CÁNCER Autor : Virginia Cabello Moreno. Director : Juan Manuel Artacho Terrer. Zaragoza, Mayo de 2012 Departamento de Ingeniería Electrónica y Comunicaciones SISTEMAS DE OPTIMIZACIÓN DIRECTA S DE MÁQUINA EN PROCESOS DE RADIOTERAPIA PARA EL TRATAMIENTO DEL CÁNCER : Juan Manuel Artacho Terrer. III RESUMEN SISTEMAS DE OPTIMIZACIÓN DIRECTA DE PARÁMETROS DE MÁQUINA EN PROCESOS DE RADIOTERAPIA PARA EL TRATAMIENTO DEL CÁNCER La Radioterapia de Intensidad Modulada (IMRT) es una técnica reciente avanzada que permite ajustar la radiación a la forma tridimensional del tumor, limitando a su vez la dosis que reciben los tejidos sanos adyacentes. Para ello se utilizan colimadores multiláminas (MLC) que administran el tratamiento sobre el paciente mediante un conjunto de aperturas o huecos a través de los cuales pasa la radiación. Por tanto un aspecto clave en IMRT es determinar las formas de las aperturas y los tiempos de exposición de cada una de ellas. Hasta el presente, el problema de la optimización se ha planteado mediante la determinación previa de las aportaciones de cada haz, modeladas a través de un mapa de intensidad o matriz de fluencia. Dicho planteamiento exige la descomposición posterior de esas matrices en aperturas y tiempos de exposición, determinando el proceso de administración de la dosis sobre el paciente. Las nuevas líneas de investigación apuntan hacia una unificación de ambas tareas, es decir, optimizar directamente las aperturas, sin tener que calcular previamente los mapas de intensidad, optimizando directamente los parámetros entregados en el tratamiento. Este proyecto aborda el planteamiento de la optimización mediante la combinación de los dos métodos mencionados anteriormente. De esta manera obtenemos las ventajas de cada uno de ellos. Por un lado la optimización se hace en una única etapa evitando los inconvenientes que supone hacerlo en dos y se impone como partida el número de etapas del tratamiento, limitando la complejidad del tratamiento. Por otro lado, al utilizar métodos de optimización condicionada, conseguimos tener un mayor control de la dosis aportada, tanto al tumor como a los órganos de riesgo. Una vez realizado el estudio tanto de las técnicas actuales como de las anteriores, se realiza el desarrollo del algoritmo en Matlab y se hace una comparación de los métodos propuestos con un proceso DAO. Los resultados obtenidos son satisfactorios, ya que con un número reducido de etapas, obtenemos planificaciones más homogéneas, en las que ninguno de los pesos de un haz determinado se apodera de los demás y las aperturas son más compactas, con los consiguientes beneficios que esto conlleva en el tratamiento entregado al paciente. IV A mis padres y hermanas por estar siempre cuando los he necesitado. A mis sobrinos por los ratos buenos que me hacen pasar y sacarme siempre una sonrisa. A mis amigas de la Universidad, por todos los ratos buenos y no tan buenos que hemos pasado juntas en todos estos años. Al grupo de trabajo, especialmente a Juanma por su ayuda y amabilidad constante. Y a ti Javi, por estar siempre ahí y apoyarme en todo lo que hago, sin ti esto no habría sido posible. Muchas gracias a todos. V ÍNDICE GENERAL RESUMEN ..................................................................................................................................... III ÍNDICE GENERAL ........................................................................................................................... V ÍNDICE DE FIGURAS..................................................................................................................... VII GROSARIO DE SIGLAS ................................................................................................................ VIII CAPÍTULO 1: INTRODUCCIÓN ...................................................................................................... 1 1.1. EL CANCER Y SUS TRATAMIENTOS .................................................................................... 1 1.2. RADIOTERAPIA DE INTENSIDAD MODULADA: IMRT ......................................................... 3 1.3. OBJETIVOS DEL PROYECTO ................................................................................................ 5 1.4. ESTADO DEL ARTE .............................................................................................................. 6 1.5. CONTENIDOS DE LA MEMORIA ......................................................................................... 6 CAPÍTULO 2: OPTIMIZACIÓN BASADA EN MAPAS DE FLUENCIA vs OPTIMIZACION DIRECTA DE APERTURAS .................................................................................................................................. 8 2.1. FASE INICIAL EN UN PROCESO DE RADIOTERAPIA ............................................................ 8 2.2. PLANTEAMIENTO DEL PROBLEMA .................................................................................. 10 2.2.1. DISCRETIZACIÓN DEL PROBLEMA ............................................................................ 10 2.2.2. MODELO DE RADIACIÓN .......................................................................................... 10 2.2.3. FUNCIÓN OBJETIVO .................................................................................................. 11 2.3. OPTIMIZACIÓN BASADA EN MAPAS DE FLUENCIA .......................................................... 13 2.3.1. PLANIFICACIÓN INVERSA ......................................................................................... 13 2.3.2. DESCOMPOSICIÓN EN APERTURAS MLC .................................................................. 14 2.4. OPTIMIZACIÓN DIRECTA DE APERTURAS ........................................................................ 16 2.4.1. ALGORITMO SIMULATED ANNEALING ..................................................................... 16 2.4.2. MÉTODO DE OPTIMIZACIÓN DE DAO ...................................................................... 17 2.4.3. CODIFICACIÓN DE LAS APERTURAS .......................................................................... 20 2.5. VENTAJAS DE LA OPTIMIZACIÓN DIRECTA DE PARÁMETROS DE MÁQUINA .................. 21 CAPÍTULO 3: MÉTODO COMBINADO FMO-DAO ...................................................................... 23 3.1. ESTRUCTURA DEL MÉTODO COMBINADO FMO-DAO ..................................................... 23 3.2. DESCRIPCIÓN DETALLADA DEL ALGORITMO ................................................................... 25 3.2.1. ELECCIÓN DE LA FUNCIÓN OBJETIVO ....................................................................... 25 3.2.2. INICIALIZACIÓN DEL ALGORITMO ............................................................................ 26 3.2.3. OPTIMIZACIÓN DE PESOS MEDIANTE MAPAS DE FLUENCIA ................................... 27 3.2.4. CONSIDERACIONES SOBRE FMO Y DAO ................................................................... 32 VI CAPÍTULO 4: RESULTADOS ........................................................................................................ 35 4.1. INTRODUCCIÓN ............................................................................................................... 35 4.2. RESULTADOS Y COMPARATIVA ....................................................................................... 38 4.2.1. RESULTADOS CON EL MODELO DE RADIACIÓN PRIMARIO ...................................... 38 4.2.2. RESULTADOS CON EL MODELO DE RADIACIÓN PENCIL BEAM ................................ 47 4.3. DISPERSIÓN DE LA FUNCIÓN OBJETIVO .......................................................................... 53 CAPÍTULO 5: CONCLUSIONES Y LÍNEAS FUTURAS .................................................................... 54 5.1. CONCLUSIONES ............................................................................................................... 54 5.2. LÍNEAS FUTURAS ............................................................................................................. 55 BIBLIOGRAFÍA Y REFERENCIAS .................................................................................................. 56 ANEXO I: MÉTODO LEMKE ......................................................................................................... 57 ANEXO II: MÉTODO ACTIVE SET ................................................................................................ 62 ANEXO III: RESULTADOS OBTENIDOS CON 5 SEGMENTOS POR HAZ ....................................... 64 VII ÍNDICE DE FIGURAS Figura 1.1: Desarrollo celular ....................................................................................................... 1 Figura 1.2: Comparación de distribución de dosis entre distintos tipos de radioterapia ............. 3 Figura 1.3: Matriz de fluencia y resultado de su descomposición en aperturas MLC. ................ 4 Figura 1.4: Acelerador lineal . ....................................................................................................... 4 Figura 1.5: Colimador multiláminas MLC . ................................................................................... 4 Figura 2.1: Histograma dosis-volumen . ....................................................................................... 9 Figura 2.2: Esquema general de radiación . ............................................................................... 10 Figura 2.3: Ejemplo de planificación con mapas de fluencia con 3 haces de radiación . ........... 13 Figura 2.4: Proceso de segmentación . ...................................................................................... 15 Figura 2.5: Restricciones del colimador multiláminas . .............................................................. 15 Figura 2.6: Ejemplo de la nueva codificación utilizada para las aperturas . .............................. 21 Figura 4.1: Posiciones angulares de los haces considerados . ................................................... 36 Figura 4.2: Proyecciones del volumen de los órganos considerados sobre el plano de haz . ... 37 Figura 4.3: Matrices de fluencia . ............................................................................................... 39 Figura 4.4: Segmentación de las matrices de fluencia . ............................................................. 41 Figura 4.5: Vista de la hoja del MLC con un diseño tongue and groove . .................................. 42 Figura 4.6: Efecto dosimétrico tongue and groove . .................................................................. 42 Figura 4.7: HDV integral . ........................................................................................................... 44 Figura 4.8: Comparativa de HDV . ............................................................................................. 44 Figura 4.9: Representación de la dosis recibida en 3D por la próstata . ................................... 45 Figura 4.10: Representación de la dosis recibida en 3D por el recto . ...................................... 46 Figura 4.11: Representación de la dosis recibida en 3D por la vejiga . ..................................... 46 Figura 4.12: Matrices de fluencia con modelo de radiación pencil beam . ............................... 48 Figura 4.13: Segmentación de matrices de fluencia con modelo de radiación pencil beam . .. 50 Figura 4.14: HDV integral con modelo de radiación pencil beam . ........................................... 51 Figura 4.15: Comparativa de HDV con modelo de radiación pencil beam . .............................. 51 Figura 4.16: Ampliación de la comparativa de HDV con modelo de radiación pencil beam . .. 52 Figura 4.17: Dispersión de la función objetivo . ........................................................................ 53 Figura I.1: Representación 3D de la función objetivo . ............................................................. 60 Figura I.2: Representación 3D de la función objetivo en la región factible. . ........................... 60 Figura I.3: Proyección sobre x 1 de la función objetivo en la región factible . ........................... 61 Figura I.4: Proyección sobre x 2 de la función objetivo en la región factible . ........................... 61 Figura III.1: Matrices de fluencia con 5 segmentos por haz . .................................................... 64 Figura III.2: Segmentación de las matrices de fluencia con 5 segmentos por haz . .................. 65 Figura III.3: HDV integral con 5 segmentos por haz . ................................................................ 66 VIII GLOSARIO DE SIGLAS SEOM Sociedad Española de Oncología Médica. IMRT Intensity Modulated Radiation Therapy (Radioterapia de Intensidad Modulada) CTV Clinical Target Volume (Volumen Blanco Clínico o Tumor) OAR Organ At Risk (Órgano en Riesgo) MLC Multileaf Colimator (Colimador Multiláminas) DAO Direct Aperture Optimization (Optimización Directa de Aperturas) DMPO Direct Machine Parameters Optimization (Optimización Directa de Parámetros de Máquina) GTC Grupo de Tecnologías de las Comunicaciones I3A Instituto Universitario de Investigación en Ingeniería de Aragón PCRT Planificación Computada para Radioterapia FMO Fluence Map Optimization (Optimización de Mapas de Fluencia) TAC Tomografía Axial Computarizada HDV Histograma Dosis-Volumen BEV Beam Eye of View (Visión Focal del Haz) LCP Linear Complementary Problem (Problema de Complementariedad Lineal) KKT Condiciones Karush-Kuhn-Tucker 1 1. INTRODUCCIÓN 1.1. EL CANCER Y SUS TRATAMIENTOS. El cáncer es un conjunto de enfermedades producidas por el crecimiento anormal de las células, las cuales pueden producir la formación de unas masas llamadas “tumores” o “neoplasias”. Cuando se forma un tumor maligno (cáncer) las células se dividen y crecen más allá de los límites normales, invadiendo tejidos sanos circundantes. En ocasiones se puede originar metástasis, la cual se produce cuando las células cancerígenas se propagan a distancia, normalmente por vía linfática o sanguínea, generando nuevos tumores en zonas apartadas. En la figura 1.1.a podemos ver el desarrollo de una célula normal, la cual se divide y muere durante un periodo de tiempo programado y en la figura 1.1.b se muestra el desarrollo de una célula cancerígena. [1] Actualmente, el cáncer es una de las enfermedades más comunes en todo el mundo y una de las principales causas de mortalidad; según datos de la SEOM (Sociedad Española de Oncología Médica), desde el 2005 los tumores ocupan la primera causa de muerte en los varones y la segunda en el caso de las mujeres. Un informe realizado por la SEOM que recopila datos tomados durante los años comprendidos entre 1975 y 2006, señala que durante el año 2006 se produjeron 61.000 muertes de varones y 37.000 de mujeres en España, de modo que, 3 de cada 1000 hombres y 2 de cada 1000 mujeres murieron en España a causa de algún tipo de cáncer. A pesar de que el cáncer se ha convertido en una de las principales causas de mortalidad, en los últimos años se ha conseguido una alta tasa de supervivencia tras el diagnostico de la enfermedad. Otro informe anterior publicado por la SEOM [2], que toma los datos del estudio Eurocare 4 (publicado en The Lancet Oncology 2007 y en European Journal of Cancer), muestra una evolución positiva de la supervivencia a causa del cáncer, tanto en hombres como mujeres españoles, siendo una tasa cercana al 50% en los hombres y del 59% en las mujeres. Figura 1.1: Desarrollo celular: (a) Célula normal, (b) Célula tumoral. 8 2. OPTIMIZACIÓN BASADA EN MAPAS DE FLUENCIA vs OPTIMIZACIÓN DIRECTA DE APERTURAS La fase de planificación de un proceso de radioterapia tiene como objetivo encontrar las aportaciones de cada uno de los haces de radicación, de manera que se cumplan las especificaciones de dosis impuestas por el radioterapeuta. Además, es necesario que el tratamiento planificado tenga un número de etapas asumible, siendo deseable que sea lo más mínimo posible, con el fin de reducir el tiempo de exposición del paciente a la radiación. Para la consecución del objetivo hay dos planteamientos fundamentales a la hora de obtener la planificación de la dosis: - Optimización de mapas de fluencia, en la que se optimizan unas matrices que representan los mapas de intensidad o fluencia, de cada uno de los haces. Posteriormente es preciso descomponer cada uno de los mapas de fluencia en un conjunto de aperturas y tiempos de exposición, que son los parámetros que configuran físicamente el MLC. - Optimización directa de aperturas, donde se optimizan directamente los parámetros del tratamiento, es decir las aperturas y sus pesos, relacionados directamente con los tiempos de exposición. Todo proceso de radioterapia consta de una parte inicial a partir de la cual se caracteriza el problema. Comenzaremos este capítulo explicando algunos conceptos claves para entender un proceso de radioterapia. 2.1. FASE INICIAL EN UN PROCESO DE RADIOTERAPIA. El primer paso en una planificación de radioterapia antes de comenzar con la fase de optimización propiamente dicha, es conocer los datos de los que partimos para resolver nuestro problema. La calidad del tratamiento de radioterapia depende de la exactitud con que se conocen los datos geométricos del paciente, para ello es necesario delimitar de la manera más precisa posible las diferentes regiones de interés: región tumoral (CTV) y estructuras críticas Capítulo 2 Optimización basada en mapas de fluencia vs optimización directa de aperturas 9 circundantes (OAR), a partir de la adquisición de las imagen médica del paciente en 3D (TAC: Tomografía Axial Computarizada). Por otro lado, el médico especialista, una vez que conoce el tipo de tumor al que nos enfrentamos, debe proporcionar las especificaciones de dosis adecuadas (D p ) que debe recibir cada uno de los órganos, tanto el tumor como los órganos de riesgo (OAR), que normalmente vienen dadas en forma de histograma dosis-volumen (HDV). Un HDV es una representación en tanto por ciento de dosis obtenida en cada órgano para el volumen que le corresponda, es decir, se representan el porcentaje de volumen de un órgano determinado que recibe X por ciento de la dosis prescrita por el especialista. El caso de un tratamiento ideal sería aquel en el que el 100% del CTV recibiera el 100% de la dosis (lo que permitiría la eliminación total del tumor) y los distintos OAR no recibieran nada de dosis en todo su volumen (no quedarían dañados o afectados). En la práctica esto es imposible de conseguir, debido a la proximidad que existe entre los distintos órganos; por lo que el especialista fija unos umbrales a partir de los cuales los tejidos sanos pueden recuperarse y toda planificación que esté por debajo de ellos será admitida. En la Figura 2.1 se puede observar un HDV tomado como punto de partida para obtener un tratamiento eficaz y válido. Se contempla la curva máxima para el OAR (curva inferior) y la curva mínima para el CTV. Las cruces representa los requisitos prescritos y las áreas sombreadas son las violaciones cometidas; en este caso no sería una planificación valida. Figura 2.1: Histograma dosis-volumen. El radioterapeuta es también el que determina el número de haces del tratamiento y las posiciones de cada uno de los focos de radiación, que se consiguen mediante la rotación del brazo del acelerador lineal y de la camilla. Capítulo 2 Optimización basada en mapas de fluencia vs optimización directa de aperturas 10 2.2. PLANTEAMIENTO DEL PROBLEMA. Una vez fijados los datos iniciales del tratamiento de radioterapia, el siguiente paso es plantear en términos matemáticos el problema de optimización o planificación, de cara a obtener las formas de las aperturas y pesos asociados a cada uno de los haces de radiación, que serán utilizados en la administración del tratamiento. 2.2.1. DISCRETIZACIÓN DEL PROBLEMA. Para poder tratar el problema matemáticamente necesitamos discretizar cada uno de los órganos en cubos llamados vóxeles. El tamaño de éstos tiene que ser tal que la representación del órgano sea de forma precisa. Además, cada uno de los haces formará una retícula donde cada cuadro se denomina bixel. El tamaño del bixel en vertical lo imponen las condiciones del MLC que estemos utilizando y en horizontal lo tendremos que determinar, teniendo en cuenta que debe ir en concordancia con el tamaño del vóxel. Figura 2.2: Esquema general de radiación. 2.2.2. MODELO DE RADIACIÓN. Una vez discretizado el problema es necesario conocer la dosis que llega a cada punto del paciente y las aportaciones de cada haz. Para ello se utiliza un sistema trazador de rayos [4] mediante el cual se conoce que vóxeles son radiados por cada foco. Todos los haces convergen en un punto al que se le denomina isocentro. Cada posición de la retícula puede radiar a más de un vóxel, lo que provoca un aumento de la complejidad. Capítulo 2 Optimización basada en mapas de fluencia vs optimización directa de aperturas 11 Para conocer el valor exacto que llega a cada vóxel se utiliza un sistema de radiación primario, el cual se modela mediante una ley exponencial; este modelo fue desarrollado en un proyecto anterior [5]. La atenuación desde el bixel i hasta el vóxel j se representa mediante   :   =      (2.1) Donde: d = distancia desde el foco hasta el isocentro h = distancia desde la piel del paciente hasta el vóxel µ = coeficiente de atenuación del medio Si consideramos todas las combinaciones posibles, podemos obtener una matriz de atenuaciones, F. La dosis que recibirá el paciente (D 0 ) se puede calcular de la siguiente manera:   =· (2.2) En particular en el vóxel i:   =∑      (2.3) Donde  es una matriz expresada en forma de vector que modelan la aportación de radiación y constituyen nuestras incógnitas a determinar. 2.2.3. FUNCIÓN OBJETIVO. Cualquier proceso de optimización requiere localizar una función objetivo que relacione la solución deseada con el problema, es decir, relacionará nuestros datos con nuestras incógnitas. El siguiente paso en un proceso de planificación de radioterapia, es encontrar lo que tenemos que optimizar. Nuestro objetivo es que la dosis recibida en el CTV sea igual o muy parecida a la dosis prescrita por el especialista (  ≈  ), que es uno de nuestros datos de partida. Por lo tanto, nuestra función objetivo podría ser aquella que minimizara el error cuadrático medio entre la dosis obtenida D i y la deseada D P . Si sólo tenemos en cuenta al CTV en la función objetivo, probablemente los OARs quedarán dañados, por lo que es interesante introducirlos en la función objetivo, considerando que   sea lo más próxima a cero (  ≈0). Si consideramos un único CTV y un único OAR, nuestra función objetivo será tal que: =  −   " +  −   " $ %&'  $ ()*  (2.4) Capítulo 2 Optimización basada en mapas de fluencia vs optimización directa de aperturas 12 Normalmente estos sumandos son multiplicados por un factor p CTV o p OAR , a los que se le denomina prioridades. Estos factores se utilizan para evitar malas planificaciones, ya que modelan la importancia relativa respecto a la optimización de unos órganos sobre otros. Para evitar la alta dependencia con estos factores se incluyen las varianzas de las distribuciones de dosis tanto del CTV como de OAR; esto es adecuado para funciones con distribuciones gaussianas, por tanto es correcto aplicarlo al CTV, pero no tanto al OAR. Además es interesante normalizar por el número de vóxeles de cada órgano, para evitar otorgándoles más importancia a aquellos órganos que son representados por un mayor número de vóxeles. Si tenemos en cuenta la ecuación 2.2, donde la dosis obtenida es igual a la matriz de atenuaciones F por una vector de pesos , podemos poner la ecuación en función de nuestras incógnitas. Además, si generalizamos la ecuación para varios CTVs y varios OARs obtenemos: =  + , N , / , $º 1   2− 3    4"+  + , N , / , $º 1   2− 3    4" $ %&'  $ ()* 5  (2.5) Nuestros datos para resolver el problema serán la matriz de atenuaciones y la dosis prescrita por el radioterapeuta y nuestras incógnitas serán la matriz de pesos W expresada en forma de vector, que representa la aportación de dosis de cada uno de los haces, en forma de matrices de fluencia. Por tanto, nos encontramos ante un problema inverso, que se podría resolver invirtiendo la matriz F, pero ésta no es cuadrada, ya que normalmente el número de vóxeles es mucho mayor que el número de pesos. Se suelen utilizar métodos iterativos para resolver este tipo de problemas, que partiendo de una situación inicial, va encontrando soluciones factibles mediante el error de la solución anterior, hasta llegar a soluciones válidas. Una vez que hemos caracterizado completamente nuestro sistema estamos en condiciones de obtener una planificación para nuestro tratamiento de radioterapia. Vamos a ver a continuación los dos planteamientos usuales de IMRT para conseguirlo. Capítulo 2 Optimización basada en mapas de fluencia vs optimización directa de aperturas 13 2.3. OPTIMIZACIÓN BASADA EN MAPAS DE FLUENCIA. La planificación basada en la determinación de mapas de fluencia se realiza en dos fases; la primera de ellas es la resolución de un problema inverso o planificación inversa y la segunda se denomina descomposición en aperturas MLC. En la figura 2.3 se muestra una imagen de lo que sería un proceso de optimización basado en mapas de fluencia con tres haces de radiación; a partir de unas matrices de fluencia se obtienen unas aperturas y una intensidad de radiación asociada a cada una de ellas, para cada uno de los haces. Figura 2.3: Ejemplo de planificación con mapas de fluencia con 3 haces de radiación. A continuación pasamos a hacer una descripción más detallada de cada una de las etapas. 2.3.1. PLANIFICACIÓN INVERSA. Una vez que hemos llegado al problema inverso que relaciona las aportaciones de cada uno de los haces con la dosis deseada, es necesario aplicar un método para resolver el problema, ya que como hemos dicho anteriormente no tiene una resolución inmediata, puesto que F no es una matriz cuadrada. Existen distintos métodos para la resolución del problema, pero nosotros nos centraremos en los algoritmos de optimización condicionada. Estos métodos permiten tener un mayor control de la dosis, mediante la imposición de restricciones impuestas tanto a los vóxeles pertenecientes a los OARs como a los pertenecientes al CTV. Los métodos de optimización condicionada serán los que posteriormente se utilicen en el método propuesto en este proyecto. Los métodos de optimización condicionada permiten imponer condiciones sobre los valores de las variables que queremos limitar. En nuestro caso dichas variables serán los vóxeles pertenecientes a los OARs, en los que limitaremos el valor de la dosis para que no sean dañados y optimizaremos lo máximo posible la dosis en el CTV. 11 12 1 2 1 2 2 2 1 1 A A A A A A A A A  A A A  A A A A M M M  w w w w w w W w w w       =         L L M M O M L 1 1 1 2 1 21 22 2 1 1 B B B B B B B B B  B B B  B B B B M M M  w w w w w w W w w w       =         L L M M O M L 11 12 1 21 22 2 1 1 C C C C C C C C C  C C C  C C C C M M M  w w w w w w W w w w       =         L L M M O M L Capítulo 2 Optimización basada en mapas de fluencia vs optimización directa de aperturas 14 Las restricciones impuestas se pueden unificar todas ellas mediante la expresión: 6∙8≤: 2.6) Donde A es la matriz de restricciones, b el vector limitación y x las incógnitas, que en nuestro caso es el vector de pesos a determinar, llamado anteriormente . Las restricciones a nuestro problema son las siguientes: - Los pesos deberán ser siempre mayores o iguales que cero, puesto que lo que representan son valores de radiación y no tiene sentido el que sean negativos. - Impondremos que los pesos estén por debajo de un cierto peso límite para evitar que unos pesos se apoderen sobre otros y produzcan quemaduras en el externo. - Además se pueden introducir restricciones para que la dosis en los OARs esté por debajo de un cierto valor, que vendrán determinada por el HDV. Es importante elegir correctamente el valor del peso límite, ya que si este tiene un valor muy elevado puede ocurrir que unos pesos se apoderen de otros, pudiendo causar quemaduras en el externo, aquellos haces que presentan mayor intensidad. Por otro lado si el peso límite es demasiado pequeño, todos los haces radiaran al valor máximo, pero será imposible conseguir las recomendaciones del especialistas para el CTV. Una vez que sabemos qué es lo que tenemos que optimizar, hemos impuesto las restricciones a nuestro problema y tenemos descrito nuestro problema inverso, queda decidir qué método de optimización condicionada será utilizado para obtener nuestros mapas de fluencia. Éstos posteriormente serán segmentados, dividiéndolos en distintas aperturas que serán las aplicadas al paciente. Entre los métodos de optimización condicionada utilizados destacaremos el método de Lemke y el método de Active Set, que serán los utilizados para la realización de este proyecto. Para obtener más información sobre estos métodos ver los anexos I y II, respectivamente. 2.3.2. DESCOMPOSICIÓN EN APERTURAS MLC. Una vez acabada la fase de planificación se obtiene como resultado un vector  del cual podemos reconstruir la matriz de fluencia, que nos permite cumplir las restricciones de dosis impuestas por el radioterapeuta. Para poder administrar esta dosis tenemos que separar la dosis que debe radiar cada haz en distintas etapas que puedan ser administradas por el colimador. Este proceso se conoce como segmentación o descomposición en aperturas MLC. A continuación se muestra un ejemplo de segmentación (Figura 2.4), donde cada matriz se corresponde con una etapa de radiación en el paciente. Capítulo 2 Optimización basada en mapas de fluencia vs optimización directa de aperturas 15 W 11 W 12 W 13 W 14 W 15 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 W 21 W 22 W 23 W 24 W 25 0 0 1 1 0 0 1 1 1 0 0 0 1 1 0 W 31 W 32 W 33 W 34 W 35 =a 1 0 0 1 1 0 +a 2 0 0 0 1 0 + …+a N 0 0 1 1 0 W 41 W 42 W 43 W 44 W 45 0 1 1 1 0 0 0 1 1 0 0 0 1 1 0 W 51 W 52 W 53 W 54 W 55 0 0 0 1 0 0 1 1 1 0 0 0 1 1 0 W 61 W 62 W 63 W 64 W 65 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 Figura 2.4: Proceso de segmentación. El proceso consiste en trasformar nuestra matriz W en una serie de segmentos o máscaras binarias, las cuales representan las aperturas del colimador. Estas máscaras son matrices de unos y ceros, donde los unos representan zonas abiertas que dejan pasar la radiación y los nulos son posiciones cerradas de las láminas donde no pasa la radiación. A demás cada matriz de unos y ceros va acompañada por un número que la multiplica y que representa la intensidad de radiación, que está directamente relacionado con el tiempo que se somete al paciente a la radiación, ya que la intensidad que genera el acelerador lineal es constante. Hay que tener en cuenta ciertas restricciones que presenta el MLC, ya que si alguna de estas se produce en las máscaras obtenidas, la solución no será válida, debido a que no será realizable físicamente. El colimador que se utiliza en el grupo de trabajo presenta las siguientes restricciones: - No podremos tener más de una apertura por fila (Figura 2.5.a). - La lámina izquierda de una fila superior o inferior no puede superar la posición de la lámina derecha de una fila y viceversa, es decir, no puede haber colisión entre una lámina izquierda y una derecha (Figura 2.5.b). - No podemos tener más de una apertura por máscara (Figura 2.5.c). Figura 2.5: Restricciones del colimador multiláminas: (a) Apertura única por fila, (b) Colisión entre láminas vecinas, (c) Múltiples aperturas por máscara. En la figura 2.5 podemos ver marcadas con un círculo rojo aperturas donde se están produciendo este tipo de restricciones. a) b) c) Capítulo 2 Optimización basada en mapas de fluencia vs optimización directa de aperturas 16 2.4. OPTIMIZACIÓN DIRECTA DE APERTURAS. La Optimización Directa de Aperturas (DAO) se basa en obtener directamente los parámetros del MLC, forma de las aperturas de las láminas de cada segmento o etapa y tiempos de exposición, que serán administradas al paciente, para conseguir la dosis deseada. La idea es obtener en un solo paso tanto las aperturas de las láminas como sus correspondientes intensidades de radiación, a partir de las cuales se obtiene la dosis prescrita mediante el HDV, que ha sido especificado por el radioterapeuta. Este planteamiento presenta ciertas ventajas, como obtener planes de radiación directamente entregables al paciente, ya que el algoritmo de optimización tiene en cuenta las restricciones que presenta el MLC y que hemos visto en la figura 2.5. De esta forma se consigue eliminar la etapa de descomposición en aperturas MLC necesaria en el caso de la optimización a través de mapas de fluencia y que en muchas ocasiones conlleva a una reducción de la calidad del tratamiento. 2.4.1. ALGORITMO SIMULATED ANNEALING. Para obtener directamente las aperturas se emplea un método meta-heurístico, llamado simulated annealing, que persigue la búsqueda de un mínimo global dentro de un problema que presenta varios mínimos locales. A diferencia de otros métodos que se quedan estancados en mínimos locales debido a que partiendo de una solución inicial sólo aceptan cambios que obtengan una mejor solución, simulated annealing acepta cambios transitorios que empeoran la función objetivo, permitiendo de esta forma no quedarse en un mínimo local. Simulated annealing se basa en el proceso de enfriamiento del acero, donde se alternan fases de enfriamiento y calentamiento para conseguir que sus cristales sean de mayores tamaños y menos defectuosos; en este estado el acero presenta un mínimo global de energía y su mayor dureza. Cuando se calienta el material los átomos se pueden mover a estados adyacentes. Se relaciona la temperatura con la probabilidad de pasar a niveles de mayor (calentamiento) o menor (enfriamiento) energía. Se parte de una temperatura alta para permitir moverse a estados de menor energía. Conforme el algoritmo va avanzando se reduce la probabilidad de moverse a un estado con menor energía, ya que la temperatura va disminuyendo paulatinamente. Esto se hace hasta que la temperatura final es cero, donde se habrá alcanzado el mínimo global. A medida que la temperatura disminuye es más difícil ir a un estado de menor energía y por lo tanto más difícil salir de un mínimo local. Capítulo 2 Optimización basada en mapas de fluencia vs optimización directa de aperturas 17 Nuestro cambio a un estado de menor energía equivale a una disminución de la función objetivo. La función objetivo utilizada en un proyecto anterior [3] del grupo de trabajo, en el que se desarrollo el método DAO, fue la que se ha presentado en la ecuación 2.5. La función objetivo nos reduce la complejidad del plan de tratamiento al minimizar un único valor numérico. La elección de la función objetivo está relacionada con la calidad del tratamiento y con los objetivos de dosis deseados. 2.4.2. MÉTODO DE OPTIMIZACIÓN DE DAO El proceso de optimización de las posiciones de las láminas del colimador MLC que configuraran cada apertura y el peso asignado a cada una de ellas, que serán nuestras variables a optimizar, parte de una inicialización de las variables utilizadas posteriormente en el algoritmo, que se repetirá hasta conseguir una solución óptima. En la fase de inicialización se les da un valor determinado a las variables que se especifican a continuación: 1. Las posiciones de las hojas de cada apertura se inicializa para radiar solamente el BEV (Beam Eye of View) de los tejidos cancerígenos. Las aperturas se adaptan a la forma que ve cada haz de los vóxeles del CTV, lo que significa que todas las aperturas del mismo haz tendrán la misma forma al comienzo del algoritmo. 2. Todos los pesos de las distintas aperturas se inicializan a un valor determinado. 3. Se selecciona una temperatura inicial. La temperatura inicial debe de ser alta, para que al comienzo sea fácil moverse a otra solución. 4. Se selecciona un esquema de enfriamiento adecuado. Los esquemas de enfriamiento más utilizados en la literatura [6] son: a. Enfriamiento exponencial: ;<= => ·;  (2.7) b. Enfriamiento rápido: ;<=  =>? ·;  (2.8) c. Enfriamiento logarítmico: ;<=  @A=>?? ·;  (2.9) Donde: · α = parámetro que escala la temperatura-tiempo. · T 0 = temperatura inicial. · n = número actual de iteraciones. El alcance o no de la solución óptima depende del esquema de enfriamiento. Considerando un esquema de enfriamiento demasiado rápido muy posiblemente llegaremos a un mínimo local, mientras que con un enfriamiento demasiado lento conseguiremos el mínimo global, pero puede provocar tiempos de búsqueda muy elevados. Capítulo 3 Método combinado FMO-DAO 24 El proceso se muestra con más detalle en el siguiente diagrama: Capítulo 3 Método combinado FMO-DAO 25 3.2. DESCRIPCIÓN DETALLADA DEL ALGORITMO. Considerando como punto de partida una optimización DAO, se proponen distintas mejoras hasta obtener el método propuesto en este trabajo. Explicamos detalladamente las partes más importantes del algoritmo utilizado y las principales aportaciones realizadas respecto al DAO original. 3.2.1. ELECCIÓN DE LA FUNCIÓN OBJETIVO. Comenzaremos eligiendo la función objetivo que será utilizada para la optimización de nuestro algoritmo. En trabajos previos del grupo de trabajo se ha trabajado principalmente con dos funciones objetivo. Una de ellas es la descrita en el apartado 2.2.3, en la que se tienen en cuenta tanto los órganos de riesgo como el tumor cancerígeno. La expresión utilizada es la siguiente: =  + , N , / , $º 1   2− 3    4"+  + , N , / , $º 1   2− 3    4" $ %&'  $ ()* 5  (3.1) Esta función objetivo presentaba dos problemas, por un lado la dependencia con las prioridades, ya que no siempre es fácil elegir la prioridad adecuada para cada órgano, tal que permita cumplir las especificaciones de dosis prescritas por el médico especialista. Por otro lado como ya hemos dicho en el punto 2.2.3, la inclusión de las varianzas en la función objetivo es más correcta para el caso del CTV que para los OARs, ya que la distribución de dosis en los OARS no responde exactamente a una gaussiana. Otra opción sería una función objetivo en la que solo aparezca el CTV. Podemos pensar que este planteamiento es erróneo porque los órganos de riesgo quedaran quemados, pero para evitar esto se introducen condiciones tales que: m  ≤m nopqrsts_ (3.2) En este caso la optimización se realiza en dos pasos, en la primera no imponemos condiciones en los OARs con el fin de detectar qué zonas de ellos recibirán más radiación. Con los resultados obtenidos podemos calcular la dosis en cada órgano y dividir el OAR en regiones según la dosis que reciban cada vóxel. Una vez realizado esto, podemos hacer un segundo proceso de optimización donde se impongan los límites de dosis a cada una de las regiones de los OAR. Capítulo 3 Método combinado FMO-DAO 26 Tras barajar las dos posibilidades y efectuar distintas pruebas con ambas, nos hemos decantado por utilizar la primera opción en nuestro trabajo, ya que, nuestro algoritmo se comporta mejor con ella. Esto sucede porque como en nuestro proceso intercalamos DAO y este utiliza la función objetivo para aceptar o no los cambios, si utilizáramos la otra opción tendríamos que estar modificando la función objetivo, ya que DAO utiliza también la primera opción. 3.2.2. INICIALIZACIÓN DEL ALGORITMO. El algoritmo comienza con una fase de inicialización de las variables que engloban tanto a las variables necesarias para el algoritmo de simulated annealing (número de iteraciones máximo, número de iteraciones permitido en cada temperatura actual, probabilidad de inicio, temperatura de enfriamiento inicial, esquema de enfriamiento, función objetivo, número de segmentos por haz), como las variables de salida (posiciones de las láminas de cada segmento y el peso asociado a cada uno de ellos). En esta fase las láminas (^ , ] ) son inicializadas con el BEV (Beam Eye of View) de los tejidos cancerígenos, o lo que es lo mismo, con los vóxeles del CTV que ve cada haz y los pesos a un valor determinado. Todos los segmentos de un mismo haz son inicializados a un mismo valor, tanto aperturas como pesos. Una vez que todas las variables tienen un valor asociado, inicializamos nuestras variables de salida (posiciones de las láminas y peso asociado de cada uno de los segmentos), que posteriormente serán optimizadas en el algoritmo iterativo, a un valor más cercano al resultado final. De esta manera mejoramos la eficiencia del método teniendo una base de partida más exacta. Se han realizado diferentes ensayos con las siguientes inicializaciones: - La primera opción fue realizar una optimización conjunta de pesos y láminas mediante DAO con unas pocas iteraciones. - La siguiente opción utilizaba un método basado en optimización de mapas de fluencia para inicializar tanto las posiciones de las láminas como sus pesos. Si elegimos la opción de inicialización mediante DAO obtenemos una mayor dosis en el CTV, al final del algoritmo. En las planificaciones obtenidas con una inicialización basada en FMO, los histogramas obtenidos para el CTV quedan muy por debajo de una planificación deseable. Por tanto, el punto de partida para nuestro algoritmo será una inicialización de los pesos y las posiciones de las láminas, mediante un proceso puramente DAO con pocas iteraciones y posteriormente antes de comenzar las iteraciones, haremos una primera optimización de los pesos mediante el método de optimización condicionada utilizado posteriormente. Capítulo 3 Método combinado FMO-DAO 27 3.2.3. OPTIMIZACIÓN DE PESOS MEDIANTE MAPAS DE FLUENCIA. A diferencia del DAO, optimizaremos los pesos mediante una técnica de optimización basada en mapas de fluencia, pero seguiremos optimizando las aperturas mediante el método utilizado en DAO. De esta forma intentamos conseguir las ventajas que presentan ambos métodos: - Por un lado un mayor control de la dosis al hacer una optimización de pesos con métodos de optimización condicionada. - Por otro lado la ventaja de DAO de poder imponer de partida el número de etapas y el hecho de no perder calidad y precisión en el tratamiento planificado. Para la elección del método de optimización condicionada utilizado en la optimización de los pesos nos hemos basado en resultados anteriores obtenidos por el grupo de trabajo, especialmente en un proyecto fin de carrera anterior [7], en el que se estudiaban distintos métodos de optimización. En dicho proyecto se llega a la conclusión de que los métodos Shor, Lemke y Active Set funcionan perfectamente con cualquier número de haces y cualquier número de vóxeles. Aunque este proyecto se decanta por el método Shor, se ha descartado su utilización ya que este método precisa ajustar otros parámetros y hacia sensibles posibles comparativas. Por tanto en este trabajo utilizaremos los métodos Lemke y Active Set para la optimización de los pesos (para ampliar información de estos métodos ver anexos I y II). Al utilizar distinto procedimiento para optimizar las aperturas y los pesos, nos surge el problema de que los dos métodos de optimización utilizan formas distintas de codificar los parámetros a optimizar, ya que uno lo hace mediante matrices que representan mapas de fluencia y el otro mediante posiciones de láminas y pesos. Por lo que precisamos poder trasladar las posiciones de las láminas en mapas de intensidad y viceversa. Para la codificación de aperturas usaremos lo descrito en el punto 2.4.3 de esta memoria. Para ello creamos un variable a la que denominaremos aperture, la cual consta de tres campos: - Segments: matriz que contiene las aperturas, mediante un vector de posiciones de dos columnas, donde cada columna representa la última posición ocupada por la lámina en cada una de las fila, según la expresión 2.13. - Weights: contiene el peso correspondiente a cada uno de los distintos segmentos. - Matrix: matriz de unos y ceros que es la traducción del campo segments a la codificación que se tenía anteriormente. Esta variable se irá actualizando en cada optimización según proceda, es decir, si solo optimizamos los pesos, se actualizara solo el campo weights, y si se optimizan las posiciones de las láminas optimizaremos tanto segments como matrix. Capítulo 3 Método combinado FMO-DAO 28 El hecho de tener dos codificaciones para representar la posición de las láminas (segments y matrix), podría parecer contraproducente para optimizar el tiempo de cómputo del algoritmo. En contrapartida, el algoritmo es mucho más rápido al utilizarlo, ya que al introducir el campo matrix podemos cambiar rápidamente de láminas y pesos a matrices de fluencia. Si multiplicamos cada matriz de unos y ceros por su peso correspondiente y sumamos el resultado de todos los segmentos, obtendremos la matriz de fluencia de cada haz, simplemente utilizando sumas y productos. De esta manera pasamos de una codificación de pesos y posiciones de láminas a matrices de fluencia, fácilmente y con poca carga computacional. Este proceso, es necesario seguirlo cuando pasamos de optimizar posiciones de láminas a optimizar pesos, pero no será necesario realizarlo al revés, ya que en nuestro proceso de optimización de pesos mediante mapas de fluencia, lo que obtenemos son directamente los pesos que se guardaran en weights. Otra de las diferencias radica en que cuando optimizamos mapas de fluencia lo hacemos de forma conjunta, es decir, se optimizan conjuntamente aperturas y pesos mediante dichas matrices, pero como hemos dicho anteriormente, en nuestro proceso de optimización solo precisamos que se optimicen los pesos. Esto hace que sea preciso reformular el problema adecuadamente. La optimización de mapas de fluencia se basa en la resolución de un problema inverso, como se explica en el apartado 2.3.1, que sigue la siguiente ecuación D = F · , donde D es la dosis deseada fijada por el médico, F es la matriz de atenuación y  son nuestras incógnitas. Ahora necesitamos que nuestras incógnitas pasen de ser un vector  a ser un vector P, que constará de los pesos que acompañan a cada máscara de cada uno de los segmentos de todos los haces. Para ello necesitamos modificar la matriz de fluencia F. A continuación veremos un ejemplo sintético, dos haces que tienen cuatro y nueve bixeles respectivamente. Cada uno de estos haces se dividirá en dos segmentos: Haz A: W 11 W 12 = a 1 0 1 +a 2 0 1 W 21 W 22 1 1 0 0 Haz B: W 11 W 12 W 13 0 1 0 0 0 1 W 21 W 22 W 23 =b 1 1 1 0 +b 2 0 1 1 W 31 W 32 W 33 1 1 0 1 1 0 Matrices obtenidas mediante DAO Capítulo 3 Método combinado FMO-DAO 29 • En un planteamiento FMO original nuestras variables serían: v=v  v " v " v "" v w v "w v xw v "w v ""w v "xw v xw v x"w v xxw   (3.3) (2 X 2) + (3 X 3) = 13 variables • Utilizando un planteamiento FMO con DAO nuestras variables serían: M= y  y " :  : "   (3.4) 2 haces X 2 segmentos por haz = 4 segmentos en total = 4 pesos = 4 variables Como se puede observar pasamos de tener 13 incógnitas a tener 4. Teniendo en cuenta que nos encontramos con un ejemplo sencillo, supone que tenemos una importante reducción en el número de incógnitas, lo cual supondrá una ventaja en el tiempo total de cómputo de nuestro algoritmo. Notar que reducimos el número de variables en el caso de la optimización de los pesos; pero no son las variables totales a determinar en nuestro problema, ya que por otro lado están las posiciones de las láminas (^ y ] ), que como ya hemos dicho se optimizan mediante un método DAO. Para poder pasar de unas incógnitas a otras, necesitamos modificar la matriz de atenuación F, de tal manera que al multiplicada por nuestro vector de pesos P, nos de cómo resultado la dosis deseada. Si seguimos el ejemplo anterior y considerando que el número de vóxeles es ocho, lo que tendríamos seria: Problema con FMO original: z { { { { |    "  x  }  ~  •  €  • ‚ ƒ ƒ ƒ ƒ „ = z { { { { { |      "   x   }   ~   •   €   •   …         "   x  "   " "  " x  " }  " ~  " •  " €  " •  " …  "   "   " "  " x  x   x "  x x  x }  x ~  x •  x €  x •  x …  x   x   x "  x x  }   } "  } x  } }  } ~  } •  } €  } •  } …  }   }   } "  } x  ~   ~ "  ~ x  ~ }  ~ ~  ~ •  ~ €  ~ •  ~ …  ~   ~   ~ "  ~ x  •   • "  • x  • }  • ~  • •  • €  • •  • …  •   •   • "  • x  €   € "  € x  € }  € ~  € •  € €  € •  € …  €   €   € "  € x  •   • "  • x  • }  • ~  • •  • €  • •  • …  •   •   • "  • x ‚ ƒ ƒ ƒ ƒ ƒ „ z { { { { { { { { { { { { | v  v " v " v "" v w v "w v xw v "w v ""w v "xw v xw v x"w v xxw ‚ ƒ ƒ ƒ ƒ ƒ ƒ ƒ ƒ ƒ ƒ ƒ ƒ „ =· (3.5) Capítulo 3 Método combinado FMO-DAO 30 Problema con FMO utilizando DAO: z { { { { |    "  x  }  ~  •  €  • ‚ ƒ ƒ ƒ ƒ „ = z { { { { { { |    †   " †   x †   } †  "  †  " " †  " x †  " } †  x  †  x " †  x x †  x } †  }  †  } " †  } x †  } } †  ~  †  ~ " †  ~ x †  ~ } †  •  †  • " †  • x †  • } †  €  †  € " †  € x †  € } †  •  †  • " †  • x †  • } † ‚ ƒ ƒ ƒ ƒ ƒ ƒ „ ‡y  y " :  : " ˆ= † ·M (3.6) Para llegar a la matriz F’ necesitamos poner el vector  en función del vector P; para ello convertimos el vector  en una suma de productos, es decir, sumamos el resultado de multiplicar la matriz de unos y ceros, obtenida mediante DAO, por su peso correspondiente, de cada uno de los segmentos: = z { { { { { { { { { { { { | v  v " v " v "" v w v "w v xw v "w v ""w v "xw v xw v x"w v xxw ‚ ƒ ƒ ƒ ƒ ƒ ƒ ƒ ƒ ƒ ƒ ƒ ƒ „ = z { { { { { { { { { { | 0 y  +y " y  y  0 :  : " :  :  +: " : " :  +: " :  +: " 0 ‚ ƒ ƒ ƒ ƒ ƒ ƒ ƒ ƒ ƒ ƒ „ = z { { { { { { { { { | 0 0 0 0 1 1 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 1 0 0 0 1 1 0 0 0 1 0 0 1 1 0 0 1 1 0 0 0 0 ‚ ƒ ƒ ƒ ƒ ƒ ƒ ƒ ƒ ƒ „ ·‡y  y " :  : " ˆ= z { { { { { { { { { | 0 0 0 0 1 1 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 1 0 0 0 1 1 0 0 0 1 0 0 1 1 0 0 1 1 0 0 0 0 ‚ ƒ ƒ ƒ ƒ ƒ ƒ ƒ ƒ ƒ „ ·M (3.7) Y entonces, de acuerdo a la información de las aperturas, obtenemos que F’ será:  † = · z { { { { { { { { { | 0 0 0 0 1 1 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 1 0 0 0 1 1 0 0 0 1 0 0 1 1 0 0 1 1 0 0 0 0 ‚ ƒ ƒ ƒ ƒ ƒ ƒ ƒ ƒ ƒ „ (3.8) Capítulo 3 Método combinado FMO-DAO 31 La nueva matriz obtenida tendrá unas dimensiones de N X (nº de segmento del conjunto de todos los haces), donde N es el número de vóxeles de nuestro problema. Teniendo en cuenta esta nueva matriz, nuestro problema de optimización inverso quedaría de la siguiente manera: = † ·‡y  y " :  : " ˆ=  † ·M (3 .9) En nuestro ejemplo sencillo vemos que el número de incógnitas pasa de 13 a 4, es decir, se reduce en más de un tercio; cuando nos encontramos en un caso real donde el número de incógnitas puede ascender hasta 10000 incógnitas, esta reducción será muy importante en la carga computacional del algoritmo. Para que el proceso de pasar de unas incógnitas a otras sea completo, necesitamos modificar la matriz de condiciones A, que contiene todas las restricciones impuestas a la función objetivo, de la misma manera que hemos modificado F. Ahora las condiciones, por ejemplo la de peso limite, ya no tiene que aplicarse a cada una de las incógnitas, sino que se deberá aplicar a una combinación lineal de ellas, es decir a la suma de alguna de ellas. Nuestra nueva matriz A’, será la matriz A multiplicada exactamente por la misma matriz de unos y ceros que hemos multiplicado F, que siguiendo con nuestro ejemplo quedaría: 6 † = 6· z { { { { { { { { { | 0 0 0 0 1 1 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 1 0 0 0 1 1 0 0 0 1 0 0 1 1 0 0 1 1 0 0 0 0 ‚ ƒ ƒ ƒ ƒ ƒ ƒ ƒ ƒ ƒ „ (3.10) Como se ha explicado anteriormente, DAO elige si se optimiza un peso o una posición de una lámina de manera estadística y posteriormente se elige que lámina se optimiza de la misma forma. Se ha optado por realizar la optimización de los pesos una vez al comienzo del proceso de optimización, justo después de la fase de inicialización y luego optimizar cada vez que se finaliza el número de iteraciones para cada temperatura de enfriamiento debido a que hacemos una optimización conjunta de todos los pesos y no de uno de ellos en particular. Capítulo 3 Método combinado FMO-DAO 32 El optimizar los pesos de esta manera nos permite eliminar una parte estadística de nuestro algoritmo, aunque no se elimina totalmente, ya que la posición de la lámina que se optimiza sigue haciéndose como en DAO y se elige estadísticamente. Parece lógico pensar que es mejor elegir lo que se quiere optimizar que elegirlo de manera estadística; existe una probabilidad mucho mayor de que la solución obtenida sea mejor si seleccionamos aquello que queremos optimizar para mejorarla. Puesto que ahora los pesos se optimizan todos de vez, se podría pensar que bastaría con hacer una optimización inicial y otra final de ellos. Hemos optado por la opción de optimizarlos al finalizar cada temperatura ya que al igual que en DAO, usamos el algoritmo simulated annealing, explicado en el punto 2.4.1, el cual acepta o no un cambio realizado si se disminuye la función objetivo, por tanto los cambios serán aceptados o no de una manera más precisa si optimizamos los pesos continuamente. Con este razonamiento, se podría pensar en optimizar los pesos cada vez que se modifica una posición de una lámina, pero esto sería computacionalmente inviable y no mejoraría mucho los resultados ya que solo se modifica la posición de una lámina y los pesos se optimizan todos a la vez. 3.2.4. CONSIDERACIONES SOBRE FMO Y DAO. 1. Número de segmentos prefijado. Los métodos de optimización de parámetros de máquina (posiciones de aperturas y pesos asociados), a diferencia de los métodos basados en FMO, tienen la ventaja de que permiten preseleccionar el número de segmentos que tendrá cada haz. Además también serán conocidos en la planificación el número de haces y sus posiciones, que son prefijados por el radioterapeuta. Prefijar el número de segmentos permite conseguir planificaciones que cumplan las condiciones que prefija el especialista, imponiendo de partida un número reducido de estos, lo que conlleva a una reducción del tiempo de exposición del paciente a la radiación y del número de etapas que se le tienen que aplicar. Esto implica un beneficio en el paciente, ya que el hecho de reducir el número de etapas permite que el tratamiento sea más eficaz, puesto que de una etapa a otra pueden influir factores no deseados, como por ejemplo que el paciente pueda moverse. Por otro lado un menor tiempo de exposición del paciente a la radiación también será favorable para su comodidad y para poder tratar un número de pacientes más elevado en un mismo periodo de tiempo. En nuestro trabajo el número de segmentos por haz se establece como un parámetro de entrada del algoritmo, aprovechando esta ventaja de los métodos de optimización DAO. El número de segmentos elegido para nuestro tratamiento, está directamente relacionado con el número de niveles de intensidad de radiación según la ecuación 2.12. Capítulo 3 Método combinado FMO-DAO 33 2. Criterio de parada y homogeneidad. El hecho de que tengamos un número de niveles de intensidad acotado, puede suponer un problema en la planificación. Recordar que el número de niveles de intensidad está directamente relacionado con el número de etapas o segmentos, según la ecuación Y=2 > − 1. En un método DAO puede suceder que las matrices de fluencia obtenidas en la planificación tengan más niveles de intensidad de los que se puedan entregar según el número de segmentos prefijado. Para evitar esto, en un proyecto anterior [3], se introdujo un criterio de parada, el cual permitía abordar la ejecución si se alcanzaba el número máximo de niveles de intensidad. Nuestro método no precisa de este criterio de parada, ya que el propio algoritmo de optimización da como resultado matrices de fluencia mucho más homogéneas, que hacen que no se supere dicho número de niveles. Este criterio también era utilizado para evitar un estancamiento de la función objetivo, es decir, que no disminuyera en sucesivas iteraciones; en nuestro caso tal estancamiento no llega a suceder, por tanto no es necesario aplicarlo por este motivo tampoco. 3. Esquema de enfriamiento. El algoritmo simulated annealing utilizado en un proceso de optimización DAO, utiliza un esquema de enfriamiento para ir de un estado a otro de menor energía, que en nuestro caso se traduce en una disminución de la función objetivo, tal y como se explica en los puntos 2.4.1 y 2.4.2 de esta memoria. El esquema de enfriamiento utilizado anteriormente era un esquema de enfriamiento de tipo logarítmico, que se corresponde con la ecuación 2.9. Al introducir los cambios explicados anteriormente al proceso DAO, nos damos cuenta que el esquema de enfriamiento tipo logarítmico no es viable para nuestro proceso, ya que los tiempos de ejecución se elevan considerablemente (por encima de doce horas). Además de otras mejoras introducidas anteriormente, como añadir un nuevo campo matrix en la variable aperture, para reducir el tiempo de ejecución, se ha optado por cambiar el esquema de enfriamiento. Utilizamos un esquema de enfriamiento rápido, que en nuestro caso seguirá la siguiente ecuación: T 0 i+1 = T 0 i +0.99 (3.11) Este esquema nos permite tener una probabilidad más alta de aceptar el cambio propuesto, así como que la probabilidad de aceptar o no el cambio, disminuya de una manera más lenta que con un esquema de enfriamiento logarítmico. Capítulo 4 Resultados 40 Además esto repercute en que las aperturas o segmentos son más compactos, como veremos en las siguientes figuras, en las cuales se puede ver la segmentación de cada una de estas matrices de fluencia. En nuestra planificación se observa que los haces de 108 O y 252 O aportan mucha más dosis que los otros tres. Esto es debido a que a través de estos dos haces existe un menor solapamiento entre CTV y OARs y permite radiar en mayor medida al tumor sin dañar excesivamente los OARs. Este hecho también se aprecia en la planificación obtenida mediante DAO. Podemos observar que los dos métodos propuestos tienen unas matrices de fluencia muy similares. A continuación, en la figura 4.4, mostramos las aperturas asociadas a las matrices de fluencia anteriores: (a) 4 5 10 15 2 4 6 8 10 12 12 5 10 15 2 4 6 8 10 12 11 5 10 15 2 4 6 8 10 12 11 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 5 5 10 15 2 4 6 8 10 12 7 5 10 15 2 4 6 8 10 12 15 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 6 5 10 15 2 4 6 8 10 12 6 5 10 15 2 4 6 8 10 12 11 5 10 15 2 4 6 8 10 12 5 5 10 15 2 4 6 8 10 12 7 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 5 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 5 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 9 5 10 15 2 4 6 8 10 12 11 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 6 5 10 15 2 4 6 8 10 12 13 5 10 15 2 4 6 8 10 12 7 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 9 5 10 15 2 4 6 8 10 12 9 5 10 15 2 4 6 8 10 12 6 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 9 5 10 15 2 4 6 8 10 12 7 5 10 15 2 4 6 8 10 12 Capítulo 4 Resultados 41 (b) (c) Figura 4.4: Segmentación de las matrices de fluencia: (a) DAO, (b) Método con optimización de pesos mediante Active Set, (c) Método con optimización de pesos mediante Lemke. Como adelantábamos al ver las matrices de fluencia obtenidas, las aperturas correspondientes a los métodos propuestos son más compactas que en el caso de DAO. El hecho de tener aperturas más compactas hace que los efectos dosimétricos no deseados sean menores. A continuación exponemos los efectos dosimétricos que hacen que la dosis realmente administrada y la planificada no coincidan exactamente: 4 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 8 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 5 5 10 15 2 4 6 8 10 12 10 5 10 15 2 4 6 8 10 12 6 5 10 15 2 4 6 8 10 12 9 5 10 15 2 4 6 8 10 12 7 5 10 15 2 4 6 8 10 12 7 5 10 15 2 4 6 8 10 12 10 5 10 15 2 4 6 8 10 12 5 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 6 5 10 15 2 4 6 8 10 12 5 5 10 15 2 4 6 8 10 12 6 5 10 15 2 4 6 8 10 12 9 5 10 15 2 4 6 8 10 12 6 5 10 15 2 4 6 8 10 12 11 5 10 15 2 4 6 8 10 12 9 5 10 15 2 4 6 8 10 12 7 5 10 15 2 4 6 8 10 12 6 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 7 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 6 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 5 5 10 15 2 4 6 8 10 12 5 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 6 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 13 5 10 15 2 4 6 8 10 12 8 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 6 5 10 15 2 4 6 8 10 12 11 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 8 5 10 15 2 4 6 8 10 12 7 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 13 5 10 15 2 4 6 8 10 12 14 5 10 15 2 4 6 8 10 12 9 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 7 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 5 5 10 15 2 4 6 8 10 12 5 5 10 15 2 4 6 8 10 12 5 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 Capítulo 4 Resultados 42 - En el MLC, para que no exista fricción entre las láminas existe una pequeña separación entre ellas, que suele ser de 0.1 mm. Esta separación hace que haya una transmisión extra de radiación entre láminas adyacentes, es decir, entre láminas superiores e inferiores. Esto puede llevar a una variación de la dosis planificada, llegando a producir una sobredosis en los órganos sanos. Este efecto es casi contra restado con un diseño tongue and groove de las láminas, aunque siempre habrá algo de fuga. En la siguiente figura se puede ver este diseño: Figura 4.5: Vista de la hoja del MLC con un diseño tongue and groove. - Este diseño evita en cierta medida la fuga inter-láminas, pero aparece el efecto tongue and groove. Se puede producir una subdosificación debida a este diseño, cuando se utilizan campos con segmentos adyacentes, por la parte lateral de las láminas. Básicamente quiere decir que la suma de las fluencias de dos segmentos no es igual a la fluencia del campo total abierto, sino que se obtiene una línea de subdosificación, ya que hay una pequeña región siempre cubierta o por el tongue o por el groove. En la figura 4.6 ponemos un ejemplo para aclarar este efecto: Figura 4.6: Efecto dosimétrico tongue and groove. - Puede haber también una transmisión extra en la zona donde se unen los extremos opuestos de las multiláminas, debido a su forma curva. En estos casos se forman unas líneas de altas dosis, entre segmentos adyacentes, conocidas como matchlines. Capítulo 4 Resultados 43 Unas aperturas irregulares hacen que existan más zonas donde se producen los efectos expuestos anteriormente, con una consiguiente variación de la dosis entregada respecto a la planificada. Esto nunca es deseable ya que se puede llegar a sobre radiar los órganos sanos o no administrar la suficiente dosis al tumor. Por otro lado, en la figura 4.4 se puede ver además de la forma de la apertura, los pesos que acompañan a cada una de ellas y que se relaciona directamente con el tiempo de exposición a la radiación. En nuestros métodos no hay ningún peso que se apodere de los demás, en un determinado haz. Además la diferencia entre el peso mayor y el menor de nuestras planificaciones, no es muy grande, a pesar de que como hemos visto anteriormente hay dos haces que aportan mayor cantidad de la radiación. Atendiendo a estos aspectos podría decirse que el método que utiliza Active Set tiene un comportamiento algo mejor que el que usa Lemke, aunque parece que las formas de las aperturas son más compactas en Lemke. En las diversas pruebas utilizadas durante este trabajo hemos observado una tendencia del DAO a que haya un segmento que se apodere de los demás dentro de un mismo haz. Este efecto no es deseado ya que se puede producir una quemadura en el externo. Una vez vistas las planificaciones obtenidas, presentamos a continuación el HDV que resulta de cada una de ellas: (a) (b) 0 20 40 60 80 100 120 0 50 100 X: 40 Y: 52.79 HDV integral % Dosis % Volumen X: 70 Y: 19.43 0 20 40 60 80 100 120 0 50 100 X: 40 Y: 57.8 % Dosis % Volumen HDV integral X: 70 Y: 21.64 Capítulo 4 Resultados 44 (c) Figura 4.7: HDV integral: (a) DAO, (b) Método con optimización de pesos mediante Active Set. (c) Método con optimización de pesos mediante Lemke. Si observamos los HDV de la figura 4.7 vemos que DAO es más conservador con los OARs, pero se observa un empeoramiento en el CTV. En un principio es mejor cuanta menos dosis se les da a los órganos de riesgo, pero siempre que se cumplan las condiciones indicadas por el radioterapeuta la planificación es igualmente válida. En nuestros métodos se cumplen las condiciones impuestas por el radioterapeuta y se aprecia una mejora en la dosis en el CTV, por tanto nuestra planificación mejoraría el HDV, aunque la diferencia no es muy notable. Para poder comparar mejor los HDV, en la siguiente figura se representan los tres HDV conjuntamente: Figura 4.8: Comparativa de HDV: Azul  Lemke, Verde  Active Set y Rojo  DAO. En la figura 4.8 las líneas azules se corresponden con una planificación que usa Lemke para optimizar los pesos, la verde es una planificación que usa Active Set y la roja representa la planificación hecha mediante DAO. 0 20 40 60 80 100 120 0 50 100 X: 40 Y: 57.7 HDV integral % Dosis % Volumen X: 70 Y: 25.42 0 20 40 60 80 100 120 0 50 100 HDV integral % Dosis % Volumen Capítulo 4 Resultados 45 Se observa que las planificaciones realizadas con nuestros dos métodos son muy similares también en cuanto a HDV se refiere. Vemos que efectivamente DAO es más conservador con los OARs, pero la dosis en el tumor es peor; además la caída de la dosis en el CTV es más abrupta en nuestro método, que también es una ventaja en la planificación, ya que hay más porcentaje del volumen del CTV que recibe más dosis. Por tanto vemos que en el HDV obtenemos una mejora con nuestro método, además de proporcionarnos matrices de fluencia más homogéneas que hace que las formas de las aperturas sean más compactas y esto favorece a los efectos dosimétricos del colimador multiláminas. A continuación podemos ver una representación en 3D de la dosis que reciben cada uno de los órganos involucrados en la planificación, la próstata (CTV), el recto y la vejiga (OARs). La dosis máxima en la planificación son 72 Gy (gray). En la parte inferior de las gráficas se puede ver una escala que representa las unidades de dosis; en la figura se representa la proporción del órgano que recibe la cantidad de dosis que apunta la flecha verde. Figura 4.9: Representación de dosis recibida en 3D por la próstata. En la figura 4.9 se observa el porcentaje del volumen de la próstata que recibe una cantidad superior a 70 Gy. Como se puede observar prácticamente toda la próstata recibe esta cantidad de dosis. Capítulo 4 Resultados 46 (a) (b) Figura 4.10: Representación de dosis recibida en 3D por el recto: (a) Para dosis próximas a 0 Gy, (b) Para una dosis de 50 Gy. En la figura 4.10 se podemos ver que prácticamente todo el recto recibe una dosis muy próxima a 0, pero sin embargo prácticamente nada del recto recibe una dosis por encima de 50 Gy. (a) (b) Figura 4.11: Representación de dosis recibida en 3D por la vejiga: (a) Para dosis próximas a 0 Gy, (b) Para una dosis de 60 Gy. Capítulo 4 Resultados 47 En la figura 4.11 vemos que ocurre algo parecido que en el caso del recto, casi toda la vejiga recibe dosis muy próximas a 0 Gy, pero prácticamente nada recibe dosis superiores a 60 Gy. Notar que en este caso la vejiga recibe cantidades de dosis más altas que el recto en un mayor volumen; esto sucede por la fisonomía de este órgano y por estar más cercano al órgano tumoral (próstata). Las dosis tridimensionales representadas anteriormente son de una planificación realizada con un método que utiliza Lemke para la optimización de pesos. No se han mostrado las gráficas pertenecientes al caso de Active Set porque al proporcionar ambos métodos soluciones tan similares, el porcentaje de dosis que recibe cada órgano es prácticamente el mismo. 4.2.2. RESULTADOS CON EL MODELO DE RADIACIÓN PENCIL BEAM. A continuación vamos a comprobar el comportamiento de los algoritmos propuestos con un modelo de radiación pencil beam, el cual aporta una descripción más correcta del comportamiento real de la radiación, ya que además de tener en cuenta la radiación sobre los vóxeles de visión directa, tiene en cuenta la dispersión de la radiación en el interior del paciente. Como hemos hecho en el apartado 4.2.1 haremos una comparativa de los dos métodos propuestos en este proyecto, optimización de pesos mediante Lemke y optimización de pesos mediante Active Set, comparándolos a su vez con el método DAO original. En la siguiente figura podemos observar los mapas de fluencia de cada uno de los tres métodos: (a) -50 0 50 -70 -65 -60 -50 0 50 -70 -65 -60 -50 0 50 -60 -55 -50 -50 0 50 -70 -65 -60 0 50 -50 0 50 -70 -65 -60 Capítulo 4 Resultados 48 (b) (c) Figura 4.12: Matrices de fluencia con modelo de radiación pencil beam: (a) DAO, (b) Método con optimización de pesos mediante Active Set, (c) Método con optimización de pesos mediante Lemke. Al igual que cuando teníamos un modelo de radiación primario, observamos que con nuestros métodos se obtienen matrices de radiación más homogéneas en todos los haces, lo cual es una ventaja, como hemos visto anteriormente, a la hora de encontrar las aperturas administradas al paciente. Los dos haces que aportaban más radiación anteriormente siguen siendo los que más aportan ahora, pero como ya hemos dicho esto es normal porque en estas posiciones existe menor solapamiento entre CTV y OARs (véase figura 4.3). -50 0 50 -70 -65 -60 -50 0 50 -70 -65 -60 -50 0 50 -60 -55 -50 -50 0 50 -70 -65 -60 0 50 -50 0 50 -70 -65 -60 -50 0 50 -70 -65 -60 -50 0 50 -70 -65 -60 -50 0 50 -60 -55 -50 -50 0 50 -70 -65 -60 0 50 -50 0 50 -70 -65 -60 Capítulo 4 Resultados 49 Por tanto vemos que aunque utilicemos un modelo de radiación más realista, nuestros métodos siguen obteniendo matrices más homogéneas que el método DAO original. Además el comportamiento es muy similar que con el modelo de radiación primario. En la siguiente figura veremos las segmentaciones obtenidas para estas matrices, de manera que podamos observar otras ventajas o inconvenientes con el modelo de radiación pencil beam: (a) (b) 22 5 10 15 2 4 6 8 10 12 5 5 10 15 2 4 6 8 10 12 8 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 5 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 5 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 21 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 10 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 4 5 10 15 2 4 6 8 10 3 5 10 15 2 4 6 8 10 3 5 10 15 2 4 6 8 10 3 5 10 15 2 4 6 8 10 3 5 10 15 2 4 6 8 10 3 5 10 15 2 4 6 8 10 9 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 9 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 14 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 6 5 10 15 2 4 6 8 10 12 12 5 10 15 2 4 6 8 10 12 6 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 11 5 10 15 2 4 6 8 10 12 7 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 6 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 5 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 6 5 10 15 2 4 6 8 10 12 9 5 10 15 2 4 6 8 10 12 6 5 10 15 2 4 6 8 10 12 9 5 10 15 2 4 6 8 10 12 7 5 10 15 2 4 6 8 10 12 7 5 10 15 2 4 6 8 10 12 5 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 5 5 10 15 2 4 6 8 10 4 5 10 15 2 4 6 8 10 3 5 10 15 2 4 6 8 10 3 5 10 15 2 4 6 8 10 4 5 10 15 2 4 6 8 10 3 5 10 15 2 4 6 8 10 11 5 10 15 2 4 6 8 10 12 9 5 10 15 2 4 6 8 10 12 6 5 10 15 2 4 6 8 10 12 8 5 10 15 2 4 6 8 10 12 5 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 6 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 6 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 5 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 56 BIBLIOGRAFÍA Y REFERENCIAS [1] Dr. E. Casado. “¿Qué es el cáncer?”, http://www.seom.org/infopublico/info-tipos- cancer/que-es-el-cancer-y-como-se-desarrolla, 2011. [2] SEOM (Sociedad Española de Oncología Médica). “El cáncer en España 2010”, http://fecma.vinagrero.es/documentos/EL_CANCER_EN_ESPA%C3%91A_2010.pdf, 2010. [3] D. Garijo. “Desarrollo de técnicas basadas en la optimización directa de aperturas aplicadas a radioterapia IMRT para el tratamiento del cáncer”, PFC de Ingeniería de Telecomunicaciones dirigido por J. Artacho Terrer, Universidad de Zaragoza, 2010. [4] J. Artacho, M. Nasarre, E. Bernués, S. Cruz. “A feasible application of constrained optimization in the IMRT system”, IEEE Trans. Biomedical Engineering, vol. 54, 370-379, 2007. [5] T. Carmona. “Ampliación de un sistema de radioterapia IMRT mediante técnicas complejas”, PFC de Ingeniería de Telecomunicaciones dirigido por E. Bernués del Río, Universidad de Zaragoza, 2003. [6] Y. Nouraniy, B. Andresenz. “A comparison of simulated annealing cooling strategies”, J. Phys. A: Math. Gen, vol. 31, 8373-8385, 1998. [7] M. Rodrigo. “Desarrollo de métodos de optimización aplicados a procesos de radioterapia de intensidad modulada”, PFC de Ingeniería de Telecomunicaciones dirigido por J. Artacho Terrer, Universidad de Zaragoza, 2008. [8] F. Carlsson. “Combining segment generation with direct step-and-shoot optimization in intensity-modulated radiation therapy”, Medical Physics, vol. 35, 3828-3839, 2008. 57 ANEXO I MÉTODO LEMKE. La aplicación de condiciones de optimización a un problema de programación cuadrática conduce a un problema lineal complementario (LCP). El método Lemke puede aplicarse para resolver problemas de este tipo. Consideramos el siguiente problema de programación cuadrática: Minimizar: " 8  ‰8+Š_‰  8+Š (I.1) Sujeto a: 68≤: (I.2) 8≥0 (I.3) Donde: ‰=2   c_Q =-2  • Š=•  • (I.4) = matriz de atenuaciones • = vector de la dosis prescrita 6 = matriz de restricciones : = vector de limitaciones Suponemos que Q es semidefinida positiva, ya que sus elementos son no negativos y, en general, alguno de ellos será nulo. La función de Lagrange correspondiente a este problema será la siguiente: •= " 8  ‰8+Š_‰  8+Š+‘  68−:+’ " +“  −8+” "  (I.5) Donde: ‘  y “  = multiplicadores de Lagrange. ’ " y ” " = variables Slack. Las condiciones KKT salen de igualar a cero la derivada de la función de Lagrange con respecto a cualquier variable que aparezca en el problema, ya sea una variable independiente, multiplicador de Lagrange o variable slack. Con lo que obtenemos las siguientes condiciones: 1) Sale de igualar a cero las derivadas parciales de la función de Lagrange con respecto a todas las variables independientes: •– •— =‰8+Š ˜ +6  ‘−“=0 (I.6) 2) Se obtiene al repetir las condiciones de ligadura del problema, o lo que es lo mismo de igualar a cero las derivadas de la función de Lagrange con respecto a los multiplicadores de Lagrange: •– •™ =−8+” " =0 (I.7) Anexo I Método Lemke 58 •– •š =68+’ " −:=0 (I.8) 3) Salen de igualar a cero las derivadas de la función de Lagrange con respecto a las variables slack: •– •› =2“  ”=0 (I.9) •– •œ =2‘  ’=0 (I.10) 4) Si aplicamos la definición de los multiplicadores de Lagrange obtendremos que, λ≥0 y u≥0. Si reescribimos la ecuación I.8 como 68+’=: se sigue cumpliendo que 2‘  ’=0. Y además de la ecuación I.7 y de la I.9 obtenemos que 2“  8=0. Si reescribimos las ecuaciones obtenemos: Optimización: ‰8+Š ˜ +6  ‘−“=0 (I.11) Factibilidad: 68+’=: (I.12) Complementariedad: 8  “+‘  ’=0 (I.13) 8≥0,’≥0,“≥0,‘≥0 (I.14) Para obtener una forma más compacta del problema agrupamos las variables de la siguiente manera: =ž“’Ÿ =ž8‘Ÿ ¡=žŠ_‰ :Ÿ Z=¢‰ 6  −6 0£ Obteniendo la siguiente descripción del problema: −Z =¡ (I.15)   =0 (I.16) Lemke introduce una variable artificial z 0 en su método, con el fin de poder construir una tabla para el sistema de ecuaciones: ¤−Z −  =¡ (I.17) Donde: ⋅ I = matriz identidad. ⋅ e = [1 1 ……. 1] T La tabla comenzara un una solución factible básica inicial, =¡+  y z=0. Con lo cual todos los multiplicadores de Lagrange y todas las variables independientes estarán igualados a cero. Dicha tabla tendrá las siguiente dimensión: (N+M) x 2(M+N+1), y tendrá el siguiente aspecto: W1 W2 W3 W4 Z1 Z2 Z3 Z4 Z0 q b ase 1=W1 base2=W2 base3=W3 base4=W4 1 0 0 0 M 11 M 12 M 13 M 14 - 1 q 1 0 1 0 0 M 21 M 22 M 23 M 24 - 1 q 2 0 0 1 0 M 31 M 32 M 33 M 34 - 1 q 3 0 0 0 1 M 41 M 42 M 43 M 44 - 1 q 4 Anexo I Método Lemke 59 Mediante un procedimiento de combinaciones lineales entre las filas y columnas de la tabla se obtiene como resultado el conjunto de variables básicas que determinan la solución. EJEMPLO: Resolvemos el siguiente problema con el método Lemke: Minimizar: 8=58  +78 " −9 " +68  +88 " −10 " Sujeto a: 38  +48 " ≤7 28  +48 " ≤5 8  ≥0,8 " ≥0 Teniendo en cuenta la ecuación I.4 y que A es la matriz de restricciones y b el vector limitación, formamos las siguientes matrices: ‰=ž122 166 166 226Ÿ Š ˜ =ž−210 −286Ÿ 6=ž3 4 2 4Ÿ :=ž75Ÿ Y por lo tanto teniendo en cuenta la definición de M y q que hemos dado anteriormente tenemos: Z=¦122 166 3 2 166 226 4 4 −3 −4 0 0 −2 −4 0 0§ ¡=¦−210 −286 75§ Como se ha indicado anteriormente la solución del problema consiste resolver la ecuación I.17. Además el algoritmo debe empezar con una solución factible inicial =¡+   y z=0, como ya habíamos indicado también. Tras varias iteraciones del método llegamos a que Z1, Z2, Z4 y W3 no serán nulas, mientras que Z3, W1, W2 y W4 sí lo serán. Podemos deducir también de la tabla final el valor de las no nulas, que serán los que aparecen en la columna q, es decir, las que están en la base: Z4=0.08 Z1=0.06 W3=1.94 Z2=1.22 Con lo que nuestra solución será: =¦  " x } §=ž8‘Ÿ=¦8  8 " ‘  ‘ " §=¦0.06 1.22 0 0.08§ =¦   "  x  } §=ž“’Ÿ=¦“  “ " ’  ’ " §=¦00 1.94 0§ La solución satisface las condiciones KKT de optimización, complementariedad y factibilidad, con lo que hemos encontrado un punto de inflexión y además es un mínimo ya que ‘≥0. Anexo I Método Lemke 60 SOLUCIÓN GRÁFICA: A continuación vamos a verificar los resultados obtenidos gráficamente. Representamos la función de error para valores de x 1 y x 2 , en escala logarítmica, es decir, representamos 10log  8  ,8 "  en lugar de 8  ,8 " . El valor mínimo se alcanza en el punto x = (-1, 2), que se corresponde con -∞ al estar representado el logaritmo. Figura I.1: Representación 3D de la función objetivo. El valor de la función de error para aquellos puntos que cumplen todas las restricciones es el siguiente: Figura I.2: Representación 3D de la función objetivo en la región factible. A continuación representamos las vistas laterales para distinguir más claramente el punto en la región factible en el que se alcanza el mínimo error. Anexo I Método Lemke 61 Figura I.3: Proyección sobre x 1 de la función objetivo en la región factible. En la figura I.3 podemos ver que el mínimo es alcanzado en x 1 =0.06. Figura I.4: Proyección sobre x 2 de la función objetivo en la región factible. En la figura anterior podemos ver que el mínimo en este caso se sitúa en x 2 =1.22. Por tanto, hemos comprobado de forma gráfica que, el mínimo condicionado se sitúa en el punto x opt = [0.06 1.22] T . Así que la solución dada por el algoritmo Lemke es correcta. 6 2 ANEXO II MÉTODO ACTIVE SET. El método Active Set resuelve problemas del tipo: Minimizar: 8= " ·8  ·‰·8+Š_‰  ·8+Š (II.1) Sujeto a: 6·8≤: (II.2) -·8=® (II.3) La matriz Q debe ser semidefinida positiva. Y por otro lado las variables que buscamos no tienen la restricción de ser mayores que cero. El algoritmo comienza con un punto inicial y el siguiente punto se calcula según la ecuación: 8 ? =8 +¯·® (II.4) Donde: ⋅ α= el tamaño del paso. ⋅ d = dirección descendente factible. En este punto puede ocurrir que alguna de las condiciones de desigualdad se cumpla como igualdad. Estas condiciones junto con las condiciones de igualdad del problema forman lo que se conoce como conjunto activo, o lo que es lo mismo Active Set. Lo primero que tenemos que hacer es encontrar la dirección descendente factible, y para ello tenemos que resolver el siguiente problema: Minimizar: " ·8 +®   ·‰·°8 +® ±+Š ˜ ·8 +®  (II.5) Sujeto a: 6Š·® =0 (II.6) Donde ⋅ Ac = matriz formada por el conjunto de condiciones activas. Si consideramos que ² =‰8 ³ +Š_‰ y reorganizamos las expresiones obtenemos que: Minimizar: " ·®  ·‰·® +²  ·® (II.7) Sujeto a: 6Š·® =0 (II.8) Estamos ante un problema de programación cuadrática con condiciones de igualdad. Las condiciones KKT para dicho problema dan lugar a la siguiente serie de ecuaciones: ‰·® +² +6Š  ·•=0 (II.9) Anexo II Método Active Set 6 3 6Š·® =0 (II.10) Donde ⋅ v = vector con los multiplicadores de Lagrange asociados a cada condición del conjunto activo. Si reescribimos ambas ecuaciones y las ponemos de forma matricial, obtenemos: ¢‰ 6Š  6Š 0£·¢® •£=¢−² 0£ (II.11) Así pues resolviendo el sistema anterior conseguiremos tanto la dirección descendente factible como los multiplicadores de Lagrange. El siguiente paso una vez conocida la dirección es encontrar el tamaño del paso. Para ello utilizaremos la siguiente fórmula: ¯=Z´µ[1. · ¸  ¸ ·—  ¸ · +y^y 6 ¹ ·®>0] (II.12) Donde: ⋅ A I y b I = matrices que contienen el conjunto de condiciones inactivas. Cuando α<1 una nueva condición de desigualdad pasa al conjunto activo en la siguiente iteración. Y por otro lado si v contiene algún elemento negativo, como son multiplicadores de Lagrange y deben ser siempre positivos, significa que alguna condición debe de extraerse del conjunto activo; si se producen varios casos se extraerá aquel que tenga un valor más negativo de v. EJEMPLO: Consideramos las siguientes matrices como datos de nuestro ejemplo: Š_‰=ž−2 1Ÿ ‰= ž1 −1 −1 2Ÿ 6=¦1 1 2 −1 −1 0 0 −1§ :=¦3400§ Tras varias iteraciones llegamos encontramos que el conjunto activo y el punto anterior, en la iteración 6, dan como resultado una dirección nulo y los multiplicadores de Lagrange son [0.111, 0.111], que son no nulos. Además en este punto se verifica en criterio de parada, con lo cual llegamos a la solución: [2.3333, 0.6667] 6 4 ANEXO III RESULTADOS OBTENIDOS CON 5 SEGMENTOS POR HAZ. En este anexo vamos a mostrar los resultados obtenidos utilizando 5 segmentos por haz. Para ello mostraremos las figuras correspondientes a las matrices de fluencia obtenidas para cada haz, las máscaras correspondientes de la descomposición MLC de dichas matrices y los HDV, para los dos métodos desarrollados en este proyecto. Los resultados que presentamos a continuación son para el caso de considerar un modelo de radiación primario. Comenzamos presentando las matrices de fluencia obtenidas tanto para el método que utiliza Lemke para la optimización de los pesos, como para el que utiliza el método de optimización condicionada Active Set: (a) (b) Figura III.1: Matrices de fluencia con 5 segmentos por haz: (a) Método con optimización de pesos mediante Lemke (b) Método con optimización de pesos mediante Active Set. -50 0 50 -60 -40 -20 0 20 40 -50 0 50 -60 -40 -20 0 20 40 -50 0 50 -60 -40 -20 0 20 40 -50 0 50 -60 -40 -20 0 20 40 0 50 -50 0 50 -60 -40 -20 0 20 40 -50 0 50 -60 -40 -20 0 20 40 -50 0 50 -60 -40 -20 0 20 40 -50 0 50 -60 -40 -20 0 20 40 -50 0 50 -60 -40 -20 0 20 40 0 50 -50 0 50 -60 -40 -20 0 20 40 Anexo III Resultados obtenidos con 5 segmentos por haz 6 5 En este caso las planificaciones obtenidas ya no son tan iguales para los dos métodos como lo eran en el caso que teníamos 7 segmentos por haz (véase figura 4.3). En este caso el método que utiliza Lemke obtiene matrices más homogéneas que el método que utiliza Active Set, sobre todo si nos fijamos en la matriz correspondiente al haz de 324 grados. (a) (b) Figura III.2: Segmentación de las matrices de fluencia con 5 segmentos por haz: (a) Método con optimización de pesos mediante Lemke (b) Método con optimización de pesos mediante Active Set. Las aperturas obtenidas son un poco menos compactas que en el caso de 7 segmentos por haz (véase figura 4.4), aunque siguen siendo más compactas que para el caso de DAO original. Las aperturas son más compactas en el caso del método que utiliza Lemke que en el caso que utiliza Active Set. En este caso sí que hay algunos segmentos que aportan más dosis que otros, pero la diferencia entre el máximo y el mínimo sigue siendo pequeña. 3 5 10 15 2 4 6 8 10 12 6 5 10 15 2 4 6 8 10 12 12 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 7 5 10 15 2 4 6 8 10 12 13 5 10 15 2 4 6 8 10 12 11 5 10 15 2 4 6 8 10 12 8 5 10 15 2 4 6 8 10 12 12 5 10 15 2 4 6 8 10 12 9 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 6 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 8 5 10 15 2 4 6 8 10 12 10 5 10 15 2 4 6 8 10 12 11 5 10 15 2 4 6 8 10 12 13 5 10 15 2 4 6 8 10 12 12 5 10 15 2 4 6 8 10 12 8 5 10 15 2 4 6 8 10 12 6 5 10 15 2 4 6 8 10 12 5 5 10 15 2 4 6 8 10 12 10 5 10 15 2 4 6 8 10 12 5 5 10 15 2 4 6 8 10 12 6 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 8 5 10 15 2 4 6 8 10 12 7 5 10 15 2 4 6 8 10 12 5 5 10 15 2 4 6 8 10 12 5 5 10 15 2 4 6 8 10 12 6 5 10 15 2 4 6 8 10 12 11 5 10 15 2 4 6 8 10 12 11 5 10 15 2 4 6 8 10 12 11 5 10 15 2 4 6 8 10 12 15 5 10 15 2 4 6 8 10 12 6 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 4 5 10 15 2 4 6 8 10 12 9 5 10 15 2 4 6 8 10 12 13 5 10 15 2 4 6 8 10 12 13 5 10 15 2 4 6 8 10 12 13 5 10 15 2 4 6 8 10 12 6 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 3 5 10 15 2 4 6 8 10 12 15 5 10 15 2 4 6 8 10 12 5 5 10 15 2 4 6 8 10 12 11 5 10 15 2 4 6 8 10 12