Fundamento Teórico - Algoritmo de Chan-Vese

Para la descripción del método de Chan-Vese, se puede asumir que una imagen I esta conformada por dos regiones de intensidades homogéneas diferentes, en este caso I ⁱ e I ^o , en donde el objeto a detectar está representado por la región I ⁱ , es decir, el área que pertenece al objeto, la cual esta delimitada por una curva C; mientras que I ^o es el área que no pertenece al objeto. El funcional de energía se puede definir tomando en consideración la región I ⁱ , como el área dentro de la curva (Cin) e I ^o como la región que esta por fuera de la curva (Cout), de la forma siguiente:

donde C es una curva variable y las constantes C1 y C2 corresponden al valor promedio de la intensidad de los píxeles de la imagen I dentro y fuera de la curva C. Cuando la curva C esta por fuera del área de interés se tiene que F1(C) > 0 y F2(C) ≈ 0, mientras que si F1(C) ≈ 0 y F2(C) > 0 son los casos en donde la curva esta dentro del área. Dicho esto, la minimización del funcional de energía ocurre cuando la curvatura esté delineando el contorno del objeto a segmentar, es decir, cuando F1(C) ≈ 0 y F2(C) ≈ 0. A este modelo se añaden algunos términos de regularización como la longitud de la curva C, y el área dentro de C, definiendo ahora el funcional de energía que se muestra a continuación:

en donde μ ≥ 0, v ≥ 0 son parámetros de regularización constantes y λ1, λ2 > 0 parámetros de peso constantes dependientes de la imagen I. Bajo este principio, la curva C se representa por medio de una curva de nivel cero de una función Lipschitz (ϕ) utilizando conjuntos de nivel, quedando la función de energía en términos de ϕ. Por otra parte, los signos que toma ésta función pueden ser representados por medio de la función Heaviside, H(ϕ):

La derivada de H(ϕ) es la función delta de Dirac (δ ₀ ), la cual se puede aproximar por:

es decir, δ ₀ es la aproximación a la función delta de Dirac, cuando ε tiende cero. Partiendo de lo anterior, el funcional de energía puede escribirse como:

Manteniendo ϕ fijo y minimizando la energía de F (C₁, C₂, ϕ) con respecto a las constantes C₁ y C₂, éstas se pueden escribir en términos de ϕ:

Por otro lado, para poder minimizar el funcional (5) en términos de la función ϕ se define su primera variación por medio de la ecuación de Euler-Lagrange. Utilizando estas consideraciones, la ecuación de diferencias se puede escribir de la forma siguiente:

Para la parte práctica e implementación del algoritmo, se utiliza la versión discreta de la ecuación (7) logrando así, la delimitación del objeto a segmentar a través de la evolución de la curva. Utilizando el descenso de gradiente, se introduce una variable artificial con respecto al tiempo de la siguiente manera:

En esta versión del método Chan-Vese se debe tomar en cuenta el número de iteraciones (n) y el tiempo entre cada una de estas (Δt). La forma discreta de la minimización del funcional de energía se muestra a continuación:

Implementación del Algoritmo de Segmentación de Chan-Vese en 3D

En este trabajo se realizaron modificaciones al modelo descrito en la sección anterior para poder realizar las segmentaciones de las estructuras de interés de forma volumétrica. Es decir, se modificaron las funciones de energía descritas anteriormente para trabajar a lo largo de los tres ejes cartesianos (x, y, z), es decir, que las funciones son generalizadas para trabajar en todo el dominio de las imágenes. La función de energía del modelo se describe como sigue:

en donde S es la superficie variable, Sin es el dominio de la región interna del objeto, Sout la región externa y Ω representa el dominio de la imagen. La primera parte de la función hace referencia a los términos de regularización considerados, en este caso el área de la superficie y el volumen interino de la superficie, los cuales fueron adecuados para trabajar con el dominio de la función de energía.

Los otros términos describen el movimiento de expansión por parte de la superficie cuando esta se encuentra adentro del área u objeto; y el movimiento de contracción cuando se encuentra por fuera del objeto. Siendo ϕ el nivel cero que representa el contorno activo, el funcional propuesto para ser minimizado es el siguiente:

Una vez derivado el funcional anterior con respecto al termino ϕ (para lograr su minimización) por medio de la ecuación de Euler-Lagrange, la ecuación de diferencias resultante se puede escribir de igual manera como la ecuación (7). Con el fin de lograr la delimitación del objeto de interés a través de la evolución de la superficie, también se utiliza la versión discreta de la ecuación de diferencias considerando el número de iteraciones y el tiempo entre cada iteración, tal y como se muestra en la ecuación (9). Para facilitar la implementación del algoritmo propuesto descrito anteriormente, en la ecuación (7) se consideró que δ ₀ (ϕ) = 1 y v = 0 basándonos en pruebas realizadas en ^[13], además de una actualización con el fin de calcular el signo del conjunto de nivel cero de la función, dando lugar a la siguiente expresión:

Los valores de los términos λ ₁ y λ ₂ de la expresión anterior, se calcularon a partir de las intensidades de los voxeles dentro y fuera de la curva. Para lograr este objetivo, se calcularon los umbrales de los pulmones y el tumor pulmonar a partir de la información de la imagen original utilizando el método de Otsu ^[7], implementado en MatLab® por medio de la función multithresh, la cual nos regresa el umbral que nos permite diferenciar los objetos o áreas de interés, en este caso los pulmones o el tumor, respecto del fondo. Una vez determinados los umbrales, se obtuvieron las varianzas de aquellos valores de la imagen que superarán dichos umbrales (fondo), como también del caso contrario, en donde los valores corresponden a alguna de las estructuras anatómicas deseadas. Así, los valores de λ's se obtienen calculando el inverso de las varianzas obtenidas.

Adicionalmente, como curva inicial en ambos casos se utilizó todo el conjunto de máscaras resultantes de una previa segmentación de las imágenes con el método de K-means, el cual es un método que se basa en el agrupamiento de clusters representados con una media ponderada de píxeles denominada centroides. El fundamento principal del algoritmo es definir K centroides (uno para cada grupo de datos de interés) dentro de la imagen a procesar ^[6]. Para esa inicialización, se utilizó la función kmeans de MatLab®, en la cual se debe especificar el número de clusters que se desean obtener. En este caso se definieron 2 clusters, lo cual permitió diferenciar la estructura anatómica de interés (pulmones o tumor pulmonar) de la demás información de la imagen (fondo y otras estructuras).

Para todas las segmentaciones realizadas, los parámetros μ y Δt tuvieron un valor de 0.1. En la segmentación de los pulmones, las constantes C₁ y C₂ fueron calculadas con el promedio de los píxeles de todo el volumen tomando en cuenta cuándo ϕ ≥ 0 y ϕ < 0 respectivamente; mientras que para el caso de la segmentación del tumor, el cálculo de las constantes se llevó a cabo con la ayuda de un subvolumen, en donde solo era posible delimitar los píxeles del objeto, facilitando así la ubicación del área en donde se encontraba la estructura.

Finalmente, el número de iteraciones utilizadas para la segmentación de los pulmones y el tumor fueron establecidas en 10 y 5 respectivamente. Estos valores fueron definidos mediante una prueba de consistencia de segmentación, en donde se buscó cual era el número mínimo de iteraciones necesarias para obtener una segmentación que defiriera en menos de 1 % respecto de una segmentación generada con un número máximo de 100 iteraciones.

Dadas las características del método propuesto, en cada iteración del algoritmo se utiliza todo el volumen de las imágenes, lo que hace que las operaciones de la función discreta del algoritmo se realicen de forma tridimensional, obteniendo de esta manera la evolución de una superficie. En resumen, la evolución del algoritmo implementado de segmentación se puede plantear utilizando la siguiente estrategia:

Gráfico 1 Estrategia para la evolución del algoritmo implementado. 

Imágenes 4DCT de Tórax

Para el análisis del algoritmo de Chan-Vese se seleccionó un conjunto de imágenes 4DCT de tórax, obtenido de la base de datos de la institución Léon Bérard Cancer Center y por el Biomedical Imaging Research Laboratory (CREATIS lab) en Lyon, Francia ^[9], en donde es posible observar estructuras pulmonares junto con sus características internas (ramificaciones bronquiales) y un tumor pulmonar de gran tamaño, localizado en el pulmón derecho. Este conjunto de datos en 4D está conformado por 10 imágenes de TAC que fueron adquiridas a un paciente en diferentes fases del ciclo respiratorio, es decir, un ciclo inspira ción-espiración es dividido en 10 imágenes de TAC. Cada una de las imágenes de TAC tienen un tamaño de 475 x 335 x 123 voxeles, con una dimensión del voxel de 0.98 mm x 0.98 mm x 2.00 mm.

Análisis Cuantitativo del Desempeño del Método de Segmentación

Para validar la precisión de las segmentaciones volumétricas obtenidas con el algoritmo implementado basado en el método de Chan-Vese, se analizaron las diferencias entre los contornos resultantes de dicho algoritmo y los contornos validados por un experto clínico. Estos contornos fueron obtenidos siguiendo la estrategia propuesta en ^[14] para el caso de la estructura de los pulmones, y de forma manual utilizando el software MIPAV ^[15] para el caso del tumor pulmonar. Las diferencias fueron cuantificadas mediante el cálculo de índices de precisión utilizados en el ámbito médico para el análisis de estructuras de interés ^[2], entre los que se pueden encontrar: el coeficiente de similitud dice (DICE), la distancia simétrica media entre contornos (DSM), y el porcentaje de distancias mayores a la dimensión del voxel ( %DMDV) ^[16].

El índice DICE es una medida de precisión global, específicamente de traslape entre las estructuras a comparar, el cual puede tener valores entre 0 (ausencia de coincidencia en las máscaras) y 1, que se traduce en una superposición perfecta entre las estructuras. Dichas estructuras están representadas por máscaras binarias generadas a partir de los contornos obtenidos en las segmentaciones. El cálculo del índice DICE se realiza de la siguiente manera:

En este caso, se calcula la intersección entre las máscaras binarias A y B (máscaras a comparar), y # representa la cardinalidad de los conjuntos.

La DSM es el valor medio del cálculo de todas las distancias Euclidianas existentes entre los puntos de dos contornos (C_A y C_B), y se estima de la siguiente manera:

donde p_A y p_B son los voxeles pertenecientes los contornos C_A y C_B, y d(p, C) es la distancia euclidiana mínima entre los voxeles y el contorno contrario a donde se encuentran. Este índice brinda una medida del error medio (desalineación) presente entre los contornos analizados.

Por último, el índice %DMDV indica el porcentaje de distancias existentes entre contornos que superan un umbral admisible de error establecido, normalmente, se utiliza la dimensión del voxel como umbral, ya que de esta manera se establece un porcentaje de las diferencias perceptibles entre las segmentaciones. Teniendo en cuenta dos conjuntos definidos de la siguiente manera:

donde DimVoxel se refiere al grosor del voxel, el índice %DMDV se puede calcular como:

Reconstrucción Tridimensional

La reconstrucción de los resultados de ambas segmentaciones (pulmones y tumor) se realizó a través de la graficación de una isosuperficie, por medio de la herramienta ya definida en el software de MatLab®.

Lo anterior consiste en calcular los valores de la superficie de los datos provenientes de un volumen, que en este caso son las máscaras resultantes de las segmentaciones, para después conectar y ajustar esos puntos en los ejes x, y, z.

En este caso, se reconstruyó cada estructura de forma individual, y además se realizó una graficación de las dos estructuras en el mismo plano, lo cual permitió observar en conjunto los dos resultados. Para lograr esto, se tuvo que graficar a modo de transparencia la estructura del exterior (pulmones), para poder visualizar el tumor interno.

El término que le da la transparencia a una isosuperficie es alpha, cuando éste término tenga un valor más cercano a 0, más traslucido se podrá observar el volumen reconstruido, mientras que en valores cercanos a 1 el volumen tomara una apariencia más opaca.

El valor de alpha utilizado para la reconstrucción de los pulmones fue de 0.1, lo que permitió visualizar no solamente el tumor pulmonar, sino además visualizar a detalle las ramificaciones internas de los pulmones.

Dinámica Pulmonar

El procedimiento de reconstrucción tridimensional se realizó en las diez fases (imágenes de TAC) del ciclo respiratorio disponibles para poder obtener información de la dinámica de las estructuras de interés a lo largo del ciclo respiratorio. En particular, para el caso del tumor pulmonar, se buscó calcular los desplazamientos máximos (en milímetros) en cada eje cartesiano.

Esto es posible de obtener utilizando las reconstrucciones del tumor pulmonar en las fases del inicio de la inspiración y el inicio de la espiración, en conjunto con la información de la dimensión del voxel (0.98 mm x 0.98 mm x 2.00 mm). Ésta información es de importancia ya que pudiera ser utilizada en la planificación de tratamientos de Tomoterapia.

Análisis de definición de iteraciones

Para determinar el número de iteraciones para las segmentaciones de los pulmones y el tumor pulmonar se realizó un análisis de robustez y consistencia del algoritmo propuesto, en el cual se comparó el resultado de una segmentación utilizando un número máximo de 100 iteraciones, con el resultado de segmentaciones usando iteraciones de 5, 10, 20, ..., hasta 100. Tomando en cuenta lo anterior, la Figura 1 muestra la diferencia de voxeles existente entre el resultado de la segmentación de los pulmones con 100 iteraciones en comparación con el resultado de las demás segmentaciones (desde 5 hasta 100). Adicionalmente, en esta Figura se puede observar una línea roja localizada en el valor de 600 voxeles, lo cual representa que existe una diferencia del 0.01 %. De esta forma, es posible apreciar que las segmentaciones son muy similares (consistentes) sin importar el número de iteraciones utilizadas en un rango entre 5 y 100 iteraciones; sin embargo, se estableció el número de iteraciones a aquel valor que tuviera un error por debajo del 0.01 %, lo cual corresponde con un número de 10 iteraciones para el caso de la segmentación de los pulmones. Para la segmentación del tumor pulmonar se realizó el mismo procedimiento para poder determinar que el número mínimo de iteraciones necesarias para una segmentación consistente fue establecido en 5.

Figura 1 Diferencia de voxeles entre segmentaciones utilizando diferente número de iteraciones. 

Es importante mencionar que aunque pareciera que el número de iteraciones no parece tener gran impacto en la precisión de la segmentación, se busca el número mínimo de iteraciones para obtener un buen resultado en el menor tiempo posible. Por lo tanto, para el caso de la segmentación tridimensional de los pulmones con 10 iteraciones se tuvo un tiempo de cómputo aproximado de 10 minutos, mientras que para el tumor pulmonar la segmentación tardó 4 minutos con 5 iteraciones. En todos los casos las segmentaciones fueron realizadas en una computadora portátil de gama media (procesador Intel Core i7 @2.5 GHz con 8 GB de memoria RAM).

Evaluación y Análisis Cuantitativo de las Segmentaciones en 3D

Para la presentación visual de los resultados obtenidos mediante el método propuesto, se utilizó la imagen de TAC en la fase final de la espiración, la cual fue considerada en este trabajo como la imagen de referencia. En la Figura 2 se puede observar un corte central de la imagen de TAC de la imagen, donde se muestran la imagen original (Figura 2 (a)) y los resultados de la segmentación de las estructuras utilizando el método propuesto; pulmones (Figura 2 (b)) y tumor pulmonar (Figura 2 (c)).

Figura 2 Resultados de la segmentación de las estructuras de interés, (a) Imagen de TAC de pulmón y tumor (fase de referencia); (b) Contorno de los pulmones obtenido por el método propuesto; (c) Contorno de tumor pulmonar obtenido por el método propuesto. 

En esta Figura es posible apreciar que el método implementado realiza una buena segmentación de ambas estructuras anatómicas de interés, ya que logra definir de manera eficiente sus contornos. En el caso de los pulmones, el método logra definir el borde externo de los pulmones, y además logra segmentar también sus características internas (ramificaciones bronquiales), que son visibles en Figura 2 (a). En cuanto al tumor pulmonar, el algoritmo propuesto logra segmenta adecuadamente la forma compleja de dicha estructura anatómica, pudiendo diferenciar pequeños detalles localizados en la periferia del tumor. Lo anterior refuerza las conclusiones reportadas en ^[8], es decir, que el método de Chan-Vese es robusto y preciso para la segmentación de estructuras anatómicas de diferentes tamaños y formas.

Para reforzar los resultados del buen desempeño cualitativo mostrado en la Figura 2, la Tabla 1 muestra los resultados provenientes del cálculo de los índices de precisión como complemento cuantitativo de la evaluación de las segmentaciones obtenidas por el método propuesto.

Tabla 1: Resultados de índices de precisión para el análisis cuantitativo de las segmentaciones de las estructuras de interés. 

En esta Tabla, se puede verificar el buen desempeño por parte del método implementado en la segmentación de las estructuras de los pulmones y el tumor. Para el índice DICE se encontraron valores cercanos a 1 (por encima de 0.90 para el caso de los pulmones y por arriba de 0.80 para el caso del tumor) lo que demuestra una gran similitud del resultado obtenido con los contornos validados por el experto clínico, sin importar la dimensión y forma de las estructuras.

En cuanto el índice de DSM, este tiene un valor de error medio de 2.55 mm para el caso de los pulmones, el cual puede deberse a los detalles (ramificaciones internas) que el método propuesto es capaz de segmentar y que en la segmentación de referencia no son considerados por el experto clínico. Lo anterior sugiere que mientras más detalles se encuentren en el interior de los pulmones, mayor será el valor de DSM y por lo tanto también el valor de %DMDV, el cual refleja el porcentaje de distancias mayores a un umbral establecido, que en este caso corresponde a 2 mm (grosor del voxel).

Sin embargo, para el caso del tumor, estos dos índices muestran valores muy cercanos a 0, lo que habla de una gran semejanza de los contornos manuales del experto con la segmentación resultante del método propuesto, a tal grado que las diferencias son menores que la dimensión del voxel.

Reconstrucción Tridimensional

En la Figura 3 se muestran las reconstrucciones tridimensionales de los resultados de las segmentaciones en la fase de referencia, en donde se puede observar de mejor forma la estructura de los pulmones (Figura 3 (a)) y del tumor pulmonar (Figura 3 (b)). Adicionalmente, se muestra la representación de ambas estructuras tridimensionales en el mismo plano (Figura 3 (c)), en donde es posible identificar el tumor de color rojo y los pulmones a modo de transparencia, lo que nos permite observar sus ramificaciones internas.

Figura 3 Reconstrucción tridimensional de las estructuras de interés, (a) Reconstrucción de los pulmones; (b) Reconstrucción en 3D del tumor pulmonar; (c) Reconstrucción de las estructuras de los pulmones y tumor en el mismo plano. 

Estas reconstrucciones nos permiten conocer a detalle las características de las estructuras de interés, como su forma y tamaño, y en el caso del tumor pulmonar, nos permite además visualizar y confirmar su ubicación dentro de los pulmones por medio de la transparencia.

Dinámica Pulmonar

Finalmente, en la Figura 4 se muestran las reconstrucciones de ambas estructuras anatómicas de interés en las diez fases del ciclo respiratorio, en donde se puede observar la fase inspiratoria (Figura 4 (a-e)), siendo 4(a) el inicio de la inspiración, y la fase espiratoria (Figura 4 (f-j)), que inicia en la 4(f).

Figura 4 Reconstrucción tridimensional de las estructuras de los pulmones y tumor pulmonar en las fases del ciclo respiratorio, (a-e) Fase inspiratoria; (f-j) Fase espiratoria. Los ejes x, y son píxeles y el eje z cortes transversales. 

Esta Figura nos permite apreciar el desplazamiento generado tanto en los pulmones como en el tumor a causa de la respiración. Para mejorar el detalle en la dinámica del tumor, en esta Figura se agregaron unas guías (líneas de color azul) que permiten comprobar el desplazamiento máximo del tumor sobre el eje z de 17 cortes transversales, que corresponden a 34.0 mm (utilizando la información de la dimensión del voxel), entre las fases de inicio de la inspiración e inicio de la espiración. Con respecto al eje x y al eje y, el tumor tuvo un desplazamiento menor, de 2.9 mm y 3.9 mm respectivamente.

Es importante mencionar que el desplazamiento máximo a lo largo del eje z era esperado, ya que de forma anatómica al movimiento que genera el diafragma en ese eje es mucho mayor que el desplazamiento de las paredes torácicas a causa de la inspiración y la espiración. Estos resultados confirman la influencia de la respiración en el movimiento de las estructuras anatómicas estudiadas, y los desplazamientos máximos obtenidos podrían ser utilizados en la planificación de tratamientos de Tomoterapia.