Obtención de imágenes elastográficas en una dimensión

UNIVERSIDAD POLITECNICA DE VALENCIA ESCUELA POLITECNICA SUPERIOR DE GANDIA I.T. Telecomunicación (Sonido e Imagen) “Obtención de imágenes elastográficas en una dimensión” TRABAJO FINAL DE CARRERA Autor: Daniel Serra Freese Directores: Francisco Camarena Femenia Francisco Javier Redondo Pastor GANDIA, 2010 2 3 Índice de contenidos Lista de tablas ............................................................................................................................................. 7 Lista de figuras ........................................................................................................................................... 9 1. Introducción .......................................................................................................................................... 15 1.1. Motivación y objetivos ................................................................................................................. 16 1.2. Organización del trabajo ............................................................................................................. 17 2. Fundamentos de los ultrasonidos ......................................................................................................... 19 2.1. Historia ........................................................................................................................................ 19 2.2. Bases físicas de los ultrasonidos ................................................................................................ 19 2.2.1. Velocidad de propagación de las ondas longitudinales .................................................... 20 2.2.2. Impedancia característica .................................................................................................. 20 2.2.3. Reflexión y refracción en una superficie plana.................................................................. 21 2.2.4. Atenuación ......................................................................................................................... 22 2.3. Elasticidad ................................................................................................................................... 22 2.3.1. Módulos y coeficientes elásticos ....................................................................................... 22 2.3.2. Deformación unitaria ......................................................................................................... 23 2.4. Ultrasonidos en pruebas no destructivas (NDT) ......................................................................... 24 2.4.1. Medición por pulsos ........................................................................................................... 24 2.5. Resolución espacial y temporal ................................................................................................. 24 2.6. Principios de formación de imagen ............................................................................................ 25 3. Principios físicos de la elastografía ...................................................................................................... 26 3.1. El tejido biológico ........................................................................................................................ 29 3.2. Comportamiento del tejido sometido a una fuerza ..................................................................... 30 3.3. Consideraciones físicas .............................................................................................................. 31 3.4. Principios de elastografía ............................................................................................................ 31 4. Procesado de elastogramas ................................................................................................................. 33 4.1. Método de block matching para imágenes ................................................................................. 34 4.2. Block matching en elastografía ................................................................................................... 34 4.3. Tamaño del kernel ...................................................................................................................... 36 4.4. Detección de la deformación....................................................................................................... 37 4.4.1. Factor de correlación ......................................................................................................... 39 4.4.1.1. Fundamentos teóricos del factor de correlación ...................................................... 39 4.4.1.2. Algoritmo de detección mediante el factor de correlación ....................................... 39 4.4.2. Función de correlación cruzada ........................................................................................ 42 4.4.2.1. Fundamentos de la función de correlación cruzada (FCC) ..................................... 42 4.4.2.1.1. Propiedad: Desplazamiento temporal de la FCC ........................................... 43 4.4.2.2. Algoritmo de detección del desplazamiento mediante la FCC ................................ 44 4.5. Determinación de la deformación, señal de deformación diferencial ......................................... 47 4.6. Otras consideraciones ................................................................................................................ 47 4 5 5. Materiales y métodos............................................................................................................................ 49 5.1. Materiales .................................................................................................................................... 49 5.1.1. Concentraciones y formas ................................................................................................. 50 5.1.2. Elaboración de la muestra ................................................................................................. 50 5.1.3. Otras consideraciones ....................................................................................................... 51 5.1.3.1. Conservante para geles ........................................................................................... 51 5.1.3.2. Dispersores .............................................................................................................. 52 5.1.3.3. Inclusiones ............................................................................................................... 52 5.1.4. Phantoms producidos ........................................................................................................ 53 5.2. Montaje experimental .................................................................................................................. 54 5.2.1. Descripción de los elementos............................................................................................ 54 5.2.2. Montaje y colocación de los elementos ............................................................................. 57 5.3. Programa de captura de datos ................................................................................................... 59 5.4. Programa de postprocesado de datos ........................................................................................ 63 5.5. Simulación numérica ................................................................................................................... 70 5.5.1. FDTD (Finite Difference Time Domain). ............................................................................ 70 5.5.2. Programa de simulación .................................................................................................... 70 6. Resultados y discusión ......................................................................................................................... 73 6.1. Creación de phantoms ................................................................................................................ 73 6.2. Montaje experimental y captura de datos ................................................................................... 75 6.2.1. Ruido eléctrico ................................................................................................................... 76 6.3. Resultados de la simulación ....................................................................................................... 76 6.4. Postprocesado de datos reales .................................................................................................. 83 6.4.1. Comparación de los resultados del programa con resultados de recuento manual ......... 83 6.4.2. Cálculo de la velocidad de propagación ............................................................................ 87 6.4.3. Elastogramas obtenidos a partir de datos reales .............................................................. 88 6.4.4. Influencia del tamaño del kernel en la determinación de los retrasos ............................ 105 6.5. Conclusiones ............................................................................................................................. 108 Bibliografía .............................................................................................................................................. 111 Anexo 1. Script simulacion.m ................................................................................................................. 115 Anexo 2. Función fdtd2d.m ..................................................................................................................... 119 6 7 Lista de tablas Tabla 3.1: Densidades, velocidades de propagación e impedancia acústica de tejidos humanos. ........... 29 Tabla 3.2: Profundidad de penetración de pulsos ultrasonidos en tejidos humanos. .................................. 29 Tabla 3.3: Atenuación de pulsos ultrasónicos en tejidos humanos ............................................................. 25 Tabla 5.1: Concentraciones de gelatinas y conservantes ............................................................................ 50 Tabla 5.2: Solución concentrada de p-hidroxybenzoato .............................................................................. 51 Tabla 5.3: Materiales probados para actuar como dispersores ................................................................... 59 Tabla 5.4: Características del transductor V392 ......................................................................................... 56 Tabla 6.1: Dimensiones y densidades del modelo simulado ....................................................................... 76 Tabla 6.2: Deformaciones (%) aplicadas para las simulaciones .................................................................. 77 Tabla 6.3: Retrasos calculados por el programa .......................................................................................... 85 Tabla 6.4: Retrasos obtenidos de forma manual ......................................................................................... 87 Tabla 6.5: Comparación de los retrasos obtenidos por el programa con los retrasos btenidos de forma manual. ......................................................................................................... 87 Tabla 6.6: Duración del recorrido acústico en muestras .............................................................................. 88 Tabla 6.7: Resumen de deformaciones, tamaños de kernel y ganancias .................................................. 89 Tabla 6.8: Resumen de los tamaños de kernel ........................................................................................... 105 8 9 Lista de figuras Figura 1.1: Ejemplo de ecografía precompresión y postcompresión de una gelatina de tres capas. .............................................................................................................................. 15 Figura 1.2: Ecografía con elastograma superpuesto. .............................................................................. 16 Figura 2.1: Ondas reflejadas y transmitidas entre dos medios. ............................................................... 21 Figura 2.2: Señal atenuada y su envolvente. ........................................................................................... 22 Figura 2.3: Funcionamiento del ecoimpulso. ............................................................................................ 25 Figura 2.4: Comparación de los pulsos emitidos de dos ecosondas ....................................................... 26 Figura 2.5: Modos de imagen. .................................................................................................................. 27 Figura 3.1: Sistema de tres muelles de la misma elasticidad .................................................................. 31 Figura 3.2: Sistema de tres muelles de diferente elasticidad. .................................................................. 32 Figura 4.1: Algoritmo básico de block matching. ...................................................................................... 33 Figura 4.2: Ejemplo de block matching en una imagen 2D. ..................................................................... 34 Figura 4.3: Ecos producidos por dos capas cercanas en un mismo phantom. ........................................ 35 Figura 4.4: Zonas validas y ambigüedades de la señal de desplazamiento acumulado. ....................... 36 Figura 4.5: Algoritmo para detectar el desplazamiento entre señales. .................................................... 40 Figura 4.6: Kernel y rango de búsqueda para una dimensión. ................................................................ 41 Figura 4.7: Representación de la deformación para el caso ideal sin ambigüedades. ............................ 41 Figura 4.8: Desplazamiento en correlograma respecto del centro. ......................................................... 43 Figura 4.9: Ejemplo para ilustrar el funcionamiento del algoritmo, resultados parciales obtenidos con el programa de postprocesado. ...................................................................... 45 Figura 4.10: Algoritmo de de detección del desplazamiento acumulado mediante la FCC. .................... 46 Figura 4.11: Ejemplo del cálculo del desplazamiento diferencial o deformación. .................................... 47 Figura 5.1: Materiales necesarios para la elaboración de phantoms. ...................................................... 50 Figura 5.2: Phantom de dos capas sometido a una pequeña deformación ............................................ 53 Figura 5.3: Phantom de cinco capas con una sometido a una gran compresión .................................... 53 Figura 5.4: De izq. a dcha.: phantom con dos capas gruesas, phantom con una capa gruesa y phantom con dos capas finas. ................................................................................. 54 Figura 5.5: Phantom con dispersores de sésamo y fibras mucilaginosas ............................................... 54 Figura 5.6: Ordenador con tarjeta NI PCI 5102. ....................................................................................... 55 Figura 5.7: 500PR vista frontal ................................................................................................................. 55 Figura 5.8: 500PR vista posterior. ............................................................................................................ 56 Figura 5.9: Transductor V392 ................................................................................................................... 56 Figura 5.10: Montaje experimental ........................................................................................................... 57 Figura 5.11: Medición de la compresión. .................................................................................................. 58 Figura 5.12: Programa de captura de datos. ........................................................................................... 59 Figura 5.13: Control de archivo. ............................................................................................................... 60 Figura 5.14: Control de recorte. ................................................................................................................ 60 Figura 5.15: Controles verticales .............................................................................................................. 60 Figura 5.16: Control de canal. .................................................................................................................. 60 Figura 5.17: Control de resolución horizontal. .......................................................................................... 61 16 Hoy en día, la elastografía está sufriendo constantes mejoras debidas a la aplicación de nuevos algoritmos para el tratamiento de la señal (p. ej. los algoritmos usados la detección de movimiento en vídeo) y la mejora en los métodos de captación se las señales RF. También han evolucionado otras técnicas elastográficas basadas en ondas electromagnéticas (rayos X) que permiten obtener resultados similares a las técnicas ultrasónicas. Una de las tendencias es poder obtener la mayor cantidad de información en tiempo real sobre las características del tejido. Figura 1.2: Ecografía con elastograma superpuesto para una muestra de tres capas. En este trabajo se presentan algunas técnicas básicas para la detección del desplazamiento temporal de las señales RF, algoritmos, experimentos realizados en el laboratorio y conclusiones. Las señales de radiofrecuencia utilizadas en este trabajo son en una dimensión, pero los métodos y conclusiones se pueden extrapolar a mapas de radiofrecuencia en dos dimensiones. 1.1 Motivación y objetivos El objetivo principal de este trabajo es obtener imágenes elastográficas unidimensionales utilizando los medios de los que se dispone actualmente en el laboratorio. Para comprender mejor las bases teóricas y poder realizar algunos experimentos se han realizado los siguientes objetivos secundarios: • Entender y repasar los principios de la ultrasonografía. • Entender los principios de la técnica elastográfica. • Aprender a utilizar la herramienta comercial de programación LabVIEW. • Creación de un programa en MATLAB que permita simular la propagación de ultrasonidos en un modelo bidimensional. • Creación de muestras de tejido (phantoms) para la realización de los experimentos. • Creación de un programa en LabVIEW para la captura de las diferentes señales. • Procesar los datos capturados con un programa en LabVIEW para obtener los elastogramas. • Realizar varios ensayos y ajustar los valors de diferentes parámetros. 17 1.2 Organización del trabajo La presente memorias Este trabajo ha sido dividida en seis capítulos. Los cuatro primeros son una introducción teórica sobre la propagación de los ultrasonidos y las técnicas de obtención de los elastogramas. El quinto capítulo trata sobre los materiales utilizados (dispositivos, programas y muestras). Finalmente en el sexto capítulo se exponen los resultados obtenidos. Capítulo 1. Introducción. Capítulo 2. Fundamentos de los ultrasonidos. Revisión de los fundamentos de los principios físicos de la ultrasonografía y de las propiedades de los materiales. Capítulo 3. Principios físicos de la elastografía. Descripción las características de los tejidos elásticos y su comportamiento cuando son sometidos a fuerzas. Bases de la elastografía y lo que representan. Capítulo 4. Procesamiento de elastogramas. Origen y explicación de los algoritmos para la detección del desplazamiento de señales. Capítulo 5. Materiales y métodos. Descripción de los materiales utilizados. Características de las muestras y su creación. Simulaciones. Montaje experimental y procedimiento de medida. Capítulo 6. Resultados y discusión. 18 19 Capítulo 2. Fundamentos de los ultrasonidos 2.1 Historia La historia de los ultrasonidos empezó con la medida de la velocidad del sonido debajo del agua. En 1826 Jean-Daniel Colladon, un físico suizo, utilizó una campana para medir la velocidad del sonido en las aguas del lago Geneva. En las últimas décadas del siglo XIX los físicos intentaron describir las bases matemáticas de las vibraciones, transmisiones y refracciones del sonido. Uno de estos físicos fue Lord Rayleigh quien escribió el libro The Theory of Sound, que se publicó en 1877 y describía las ondas acústicas mediante ecuación matemática, creando así la base de la acústica para trabajos posteriores. Anteriormente, el biólogo italiano Lazzaro Spallanzani, demostró en 1794 que los murciélagos volaban sin problemas por la oscuridad gracias a los ultrasonidos. Más tarde el científico Francis Galton inventó el silbato que emitía ultrasonidos. La gran revolución de los ultrasonidos se produjo con el descubrimiento del efecto piezo-eléctrico en algunos cristales descubiertos por Pierre Curie y su hermano Jacques Curie en 1880. Durante la Segunda Guerra Mundial se produjeron grandes desarrollos en los sistemas de RADAR (Radio Detection and Ranging). Estos radares eran los precursores de las aplicaciones médicas en 2 dimensiones que aparecieron a finales de la década de los 40. En la década de los 50 los ultrasonidos se utilizaban más en el ámbito del tratamiento que no en la detección. Se utilizaron estos dispositivos para generar calor en los tejidos en el ámbito de la medicina de rehabilitación, p. ej. tratamiento de artritis reumática. En 1957 se desarrolló el Pan-Scanner, un aparato que rotaba alrededor del paciente sumergido parcialmente en una solución salina. Estas fueron las primeras imágenes en modo B que se obtuvieron. En las décadas de los 60 a los 70 se empezaron a producir las primeras imagines con escalas de intensidad además de los transductores transvaginales y transrectales. En esa época apenas se realizan exploraciones ultrasónicas por vía transcutánea. 2.2 Bases físicas de los ultrasonidos En un medio elástico se pueden propagar ondas mecánicas a varias frecuencias. Se llaman ultrasonidos aquellas vibraciones mecánicas que están por encima de la banda audible (20.000 Hercios). Definimos como período el tiempo que tarda un cierto patrón en repetirse, le designaremos con la letra T y su unidad fundamental son los segundos (s). La frecuencia, designada por la letra f, es la inversa del período y se mide en hercios. 1 f T = (2.1) 20 Si el sonido se propaga por un medio con una velocidad c, entonces la longitud de onda se relaciona con el período y la frecuencia con esta fórmula c cT T λ = = (2.2) Con la transmisión de la onda ultrasónica por el medio, se crea una presión acústica. Esta presión acústica está relacionada con la velocidad de las partículas del medio, designada con la letra v, mediante la siguiente fórmula p Zv = (2.3) donde Z es la impedancia acústica específica o característica (ver punto 2.2.2). 2.2.1. Velocidad de propagación de las ondas longitudinales En los sólidos se pueden propagar varios tipos de ondas con diferentes velocidades. En un medio ilimitado, p. ej. un medio en el que las dimensiones transversales son mucho mayores que la longitud de onda, la velocidad de propagación viene dada por la siguiente fórmula -1 [m s ] p k E c ρ = ⋅ (2.4) donde E p es la contante elástica y ρ es la densidad del medio. 2.2.2 Impedancia característica Como en electricidad, se puede establecer una impedancia característica del medio. Así el término c ρ es la impedancia característica del medio, puede ser real o complejo. Z c R jX ρ = = + (2.5) donde ρ es la densidad del medio y c la velocidad de propagación. Sus unidades son N·s·m -3 . Cuando una onda acústica atraviese dos medios con impedancias acústicas diferentes, se formarán ecos. Según sea la magnitud del cambio de impedancia la amplitud de la onda reflejada será mayor o menor. 21 2.2.3 Reflexión y refracción en una superficie plana La propagación de los ultrasonidos está influenciada por cualquier superficie de contorno. En esta superficie de contorno el eco será reflejado en mayor o menor medida. Si la superficie de incidencia es mucho mayor que la longitud de la onda acústica, entonces es posible aplicar las leyes de reflexión y refracción que se aplican en óptica. El factor de reflexión viene dado por la siguiente fórmula 2 02 01 0 02 01 ref inc IZ Z RI Z Z   − = =   +   (2.6) donde I inc es la intensidad de la onda incidente y I ref la intensidad de la onda reflejada. Ambas intensidades se relacionan con la intensidad transmitida (I trans ) mediante la siguiente relación ref inc trans I I I = − (2.7) En cuanto a la incidencia oblicua de ondas ultrasónicas en dos medios con diferente impedancia, el ángulo de la onda reflejada es igual al ángulo de la onda incidente y la onda transmitida será refractada según la relación entre las velocidades de ambos medios, igual como sucede con la ley de Snell. 1 2 sin sin c c α β = (2.8) donde c 1 y c 2 son las velocidades de propagación en los diferentes medios. Un caso particular es cuando se produce reflexión total, esto sucede cuando: 1 2 sin 1 c c α > (2.9) Z 2 Z 1 α β Onda incidente Onda reflejada Onda transmitida Figura 2.1 : Ondas reflejadas y transmitidas entre dos medios. 22 2.2.4 Atenuación La propagación de ondas ultrasónicas se puede describir como las deformaciones elásticas de las diferentes partículas de un medio. Dada la fricción y la conductividad térmica, estas deformaciones se ven acompañadas por pérdidas de energía de oscilación, la cual se convierte en energía calorífica. La atenuación en ultrasonidos es la reducción de la amplitud de la señal ultrasónica. Este efecto puede ser importante ya que puede afectar a la calidad de la imagen. Si se conoce la atenuación del medio se pueden introducir factores de corrección para compensar los efectos. [ ] ( ) [ ] [ ] Atenuacion dB dB MHz cm l cm f MHz α =⋅⋅ ⋅   (2.10) Las unidades de la atenuación normalmente no son lineales y se pueden tomar como referencia algunos valores representados en la tabla 3.3. Figura 2.2: Señal atenuada y su envolvente. 2.3 Elasticidad Uno de los temas más importantes de este trabajo es la elasticidad, con lo que conviene tener bien claro los conceptos más importantes que definen los medios elásticos y que nos permitirán establecer simplificaciones de los casos. En física e ingeniería, el término elasticidad designa la propiedad mecánica de ciertos materiales de sufrir deformaciones reversibles cuando se encuentran sujetos a la acción de fuerzas exteriores y de recuperar la forma original si estas fuerzas exteriores se eliminan. 2.3.1 Módulos y coeficientes elásticos. Las constantes elásticas que reciben el nombre de módulo elástico son las siguientes: Módulo de Young. Se le designa por la letra E. Está asociado directamente a los cambios de longitud que experimenta un cuerpo cuando es sometido a una fuerza de compresión o tracción. También se le conoce como módulo elástico longitudinal. 23 Módulo de compresibilidad. Se le designa por la letra K. Está asociado a los cambios de volumen que experimenta un cuerpo bajo la acción de las fuerzas (generalmente compresoras) que actúan perpendicularmente a sus superficies. No implica cambio de forma, solamente de volumen. Módulo elástico transversal. Se le designa por G. Está asociado con el cambio de forma que experimenta un material bajo la acción de fuerzas cortantes. No implica cambios de volumen, tan solo de forma. También se le conoce como módulo de cizalla, módulo elástico tangencial o módulo elástico cortante. Coeficiente de Poisson. Es coeficiente del alargamiento longitudinal entre el acortamiento de una longitud situada en un plano perpendicular a la dirección de la fuerza aplicada. Este valor coincide igualmente con el cociente de deformaciones, de hecho la fórmula usual para el coeficiente de Poisson es: trans long v ε ε = − (2.11) Para un material isótropo elástico perfectamente incompresible, este es igual a 0.5. La mayor parte de los materiales prácticos en la ingeniería rondan entre 0.0 y 0.5. En el Sistema Internacional de Unidades, los módulos se expresan en Newtons partido metro cuadrado (N/m 2 ) y el coeficiente es adimensional. 2.3.2 Deformación unitaria La magnitud más simple para medir la deformación es lo que en ingeniería se llama deformación axial o deformación unitaria, se define como el cambio de longitud por unidad de longitud 0 0 ' l l l l l ε ∆ − = = (2.12) donde l 0 es la longitud inicial de la zona en estudio y l' la longitud final o deformada. Es útil para expresar los cambios de longitud de un cuerpo o capa. 24 Cuando la deformación puede tener lugar según diversos modos y en diversas direcciones, y puede además provocar distorsiones en la forma del cuerpo, en esas condiciones se dice que el cuerpo es anísotrópico. La deformación de un cuerpo se puede caracterizar por un tensor (más exactamente un campo tensorial) de la forma [ ] 11 12 13 21 22 23 31 32 33 D ε ε ε ε ε ε ε ε ε     =       (2.13) Cuando la deformación se produce por igual en todas las direcciones se dice que el cuerpo es isotrópico. 2.4 Ultrasonidos en pruebas no destructivas (NDT) Si nos referimos al ámbito industrial, podríamos definir las pruebas no destructivas (Non-Destructive Testing) como aquellas pruebas basadas en la detección de cambios causados por irregularidades en el material (fisuras, roturas, cavidades, agujeros, inclusiones,…) bien sea por campos electromagnéticos, cambios visuales y como no, por pruebas ultrasónicas. Los métodos de prueba que utilizan ultrasonidos tienen numerosas ventajas dada su capacidad para determinar la posición y si es necesario su forma (ultrasonografía clínica). Otra de las ventajas no menos importante, es que la velocidad de propagación de los ultrasonidos es mucho menor que la velocidad de las ondas electromagnéticas, con lo cual su longitud de onda es mucho menor. Es en los tejidos heterogéneos donde a los ultrasonidos presentan una desventaja, y es que su velocidad de propagación depende enormemente de las condiciones del medio. Además, en algunos casos pueden sufrir una fuerte atenuación a lo largo de su propagación, como es el caso de la propagación en el aire. En algunos casos se requiere una mayor profundidad de penetración en el material o una mejor resolución espacial. Es allí cuando se hacen necesarias nuevas técnicas de análisis de la señal, como puede ser la manipulación del espectro o criterios de correlación. Los métodos para realizar las pruebas se dividen fundamentalmente en dos grupos. Por un lado aquellos métodos que utilizan una señal continua (armónica, periódica) y por otro lado, las que hemos utilizado para este trabajo, señales discontinuas (pulsos). 2.4.1 Medición por pulsos Este método es el que da la mayor cantidad de información sobre el material examinado, tanto desde el punto de vista de detección de irregularidades como de las dimensiones estructurales. Para que este método funcione, será fundamental la presencia de estructuras que produzcan reflexiones de las ondas ultrasónicas. Estos materiales son llamados ecogénicos. 25 En el caso de las medidas por pulsos ultrasónicos, se evalúa la duración de la propagación de la onda longitudinal ultrasónica desde el transmisor (generador de pulso) hasta el receptor (detector de pulso). El tiempo de propagación depende de la distancia entre el emisor, el objeto y el receptor; además de la velocidad de propagación y de las características del medio. En el caso de la detección por pulsos podemos decir que es análogo a la radiolocalización. En algunos documentos se refiere a este método como ecolocalización. Para detectar a qué distancia se encuentra un objeto, se considera el tiempo que el pulso tarda en llegar del emisor a un objeto y el tiempo en que la reflexión tarda en llegar al receptor. Suponiendo que el emisor del pulso y el receptor se encuentran en un mismo punto bastará con dividir el tiempo total del recorrido entre 2. 2 t d c ∆ = (2.14) El principio de funcionamiento se describe en la figura 2.3. El transductor se excita con un generador de pulsos. El receptor acoplado acústicamente al material detecta los pulsos reflejados por las irregularidades. Figura 2.3: Funcionamiento de la técnica de ecoimpulso. 2.5 Resolución espacial y temporal La resolución espacial se refiere a la finura de detalles visibles en una imagen: cuanto menor es el área representada por cada píxel en una imagen digital mayores son los detalles que pueden ser captados y mayor es la resolución espacial. En el caso de los ultrasonidos, la resolución espacial vendrá dada por el ancho del pulso y la frecuencia. Desde un punto de vista eléctrico, se puede generar un pulso con buenas características (tiempo de subida, amplitud, anchura…), pero a la hora de ser transmitido por el transductor este pulso cercano al ideal será distorsionado. Dadas las características del transductor tendremos una frecuencia de resonancia a la que oscilará durante un tiempo el transductor y una amortiguación que determinará el ancho final del pulso emitido. La frecuencia de resonancia y la velocidad longitudinal de propagación del pulso ultrasónico determinarán la longitud de onda de la señal. 32 En el siguiente caso se aumentará la constante elástica o dureza de uno de los muelles haciéndolo casi incompresible, en este caso el muelle central. Al volver a aplicar la misma fuerza que en el caso anterior podremos observar que hora si aumentamos la constante elástica de un muelle y volvemos a aplicar la fuerza, la deformación ya no será igual para todos los muelles. Los muelles con menor constante elástica se habrán deformado más que el muelle con la mayor constante elástica. Tras aplicar la fuerza y obtener el diagrama de deformación observaremos que la deformación en la zona comprimida es cero o casi cero. Esta información de deformación nos dará la información necesaria para realizar el elastograma. Este segundo ejemplo representa un material elástico, heterogéneo y isotrópico. Aunque este ejemplo sea muy sencillo e ideal, describe muy bien lo que ocurre en el mundo real. En el ámbito medicinal y para aplicaciones más avanzadas será necesario obtener mapas de deformación a partir de imágenes bidimensionales. Figura 3.2: Sistema de tres muelles de diferente elasticidad. Para aplicaciones más avanzadas será necesario tener en cuenta varios factores y correcciones como pueden ser la corrección del campo de deformaciones o el desplazamiento lateral del material sometido a una fuerza. 33 Capítulo 4. Procesado de elastogramas La mayoría de los algoritmos para el procesado de elastogramas se utilizan en aplicaciones de tratamiento digital de imágenes. Estos algoritmos se utilizan mucho en las aplicaciones de compresión de video para detectar el desplazamiento de bloques y no tener que redefinirlos otra vez. Otro ámbito de aplicación del block matching es en la detección de patrones de señal, p. ej. señales de radar. En este capítulo explicaremos la técnica utilizada para determinar el desplazamiento de las señales RF y representar los mapas de deformación. 4.1 Método de block matching para imágenes En nuestro caso, esta técnica ha sido adaptada para detectar el movimiento en un tejido cuando es sometido a una pequeña presión. La técnica de block matching se basa en algoritmos de búsqueda, comparando un patrón con otro patrón y estableciendo si se corresponden o no. En este caso los vectores de movimientos indican cuanto se ha desplazado un bloque (kernel) de nxm de una posición i a otra posición j de otra señal. Una vez obtenidos los vectores de desplazamiento ya se realizará con ellos las operaciones oportunas. Normalmente se establece una región de interés (ROI, Region of Interest) para evitar ambigüedades entre dos patrones y acelerar la búsqueda. Existen numerosas funciones para calcular el desplazamiento: • Función de correlación cruzada normalizada (NCC, Normalized Cross Correlation). • Suma de diferencias al cuadrado (SSD, Sum of Squared Differences). • Suma de diferencias absolutas (SAD, Sum of Absolute Differences). • Pel Difference Classification (PDC) Obtener frame k+1 Obtener frame k Divide en kernels de nxm Divide en kernels de nxm Obtener el desplazamiento entre kernels Vector desplazamiento Figura 4.1 : Algoritmo básico de blockmatching. 34 En la siguiente figura se puede observar el proceso de identificación de patrones mediante la técnica de block matching. La imagen de entrada k es dividida en varias regiones del cual se toma un bloque o kernel. El siguiente paso consiste en identificar la nueva posición del kernel mediante alguna de las funciones anteriores. Finalmente se obtiene la nueva posición del bloque y con ello el vector de desplazamiento. Figura 4.2: Ejemplo de block matching en una imagen 2D. 4.2 Block matching en elastografía Para determinar el desplazamiento de los componentes de dos mapas de valores será necesario al menos de un patrón y de una zona donde buscar (ver punto anterior). Para nuestro caso, se hace necesario que obtengamos como mínimo dos señales, una antes de producirse el desplazamiento o deformación y otra después de que este se produzca. La detección del desplazamiento será más facil cuanto más se parezcan las señales entre sí. Si las dos señales han variado (desplazado) mucho los resultados devueltos por los métodos de block matching pueden ser erróneos y producir ambigüedades. Esto es una limitación que se impone a la hora de capturar 35 la señal post deformación, normalmente se aconseja que este valor no supere el 5% [3]. A mayor deformación, mayores serán los errores de localización. Como se ha comentado anteriormente, las señales se dividirán en bloques o kernels (núcleos, ventanas o patrones) y se realizará la comparación entre los dos bloques. Si la forma de onda en un kernel es única, esta comparación ofrece una medida de desplazamiento igualmente única. Si la región de interés o de búsqueda es demasiado extensa puede que existan patrones más parecidos que el que realmente debería haber sido detectado. Para este caso, en nuestro trabajo hemos descartado posibles ambigüedades mediante un enventanado de Hanning. En el capítulo de discusión de este trabajo, se exponen los resultados obtenidos para los diferentes tamaños de kernels. Dado que en nuestro caso disponemos una señal dependiente de una sola variable f(t) extrapolaremos el método de block matching a una dimensión. La posibilidad de que se produzcan ambigüedades en una dimensión es mucho mayor que en dos, ya que la unicidad de los conjuntos es menor. Pongamos el caso de un pulso ultrasónico que se propaga a través de un material, si este se encuentra con un cambio de impedancia devolverá un eco. Si ahora buscamos en la señal de ecos el pulso reflejado, lo encontraremos fácilmente. Si en lugar de un único eco se forman dos y además están muy cercanos, puede que se identifique el segundo eco como el primero, dando como resultado una información errónea. Suponiendo que estamos ante una señal de postdeformación, con una deformación muy pequeña (como máximo 5%), si tomamos una muestra de la señal predeformación y buscamos el mismo patrón en la señal postdeformación, sabemos que bastará con buscar en las cercanías. A esta zona la llamaremos Región de interés (ROI, Region of Interest). Esta suposición ya nos reduce de manera importante el número de operaciones a realizar y los errores producidos por coincidencias con zonas muy lejanas que no tienen ninguna relación con el patrón pero que se parecen. En el caso de las señales obtenidas en una dimensión, bastará con algunas funciones matemáticas que nos definan la correspondencia y el desplazamiento de la señal (ver sección 4.1). Figura 4.3: Ecos producidos por dos capas cercanas en un mismo phantom. En la bibliografía consultada [3] mencionan métodos para reducir el coste computacional y obtener imágenes a tiempo real. Algunas de ellas consisten en múltiples resoluciones o en la interpolación de los 36 resultados. En nuestro caso, el coste computacional no es importante ya que se trata de una aplicación que no trabaja en tiempo real y con una cantidad muy pequeña de información. 4.3 Tamaño del kernel Escoger el tamaño adecuado del kernel no es una tarea trivial. En general kernels grandes son menos sensibles a ruido, mientras que kernels pequeños ofrecen mejores contornos entre estructuras con diferentes elasticidades. Además, cuantos mayores sean el kernel y el rango de búsqueda, también será mayor el coste computacional. Uno de los problemas más comunes es cuando surgen ambigüedades producidas por mismas estructuras acústicas. Un ejemplo para ilustrar el uso del tamaño del kernel es en la compresión de videos. En algunos casos cuando hay poco contraste entre zonas, a la hora de comprimir, las zonas se codifican y erróneamente se les atribuye un falso vector de movimiento. El resultado es una imagen de video donde se pueden ver zonas que no se mueven aun habiendo movimiento. No se debe confundir esto con el término de pixelar. Algunas consideraciones sobre el tamaño de los kernels: • Deberemos escoger un tamaño de kernel que se ajuste al tamaño de los sucesos o patrones identificables. • Deberemos escoger entre un tamaño grande que produzca pocas ambigüedades y un tamaño pequeño que ofrezca buena resolución. • Un kernel no debe abarcar más de dos sucesos similares (kernel demasiado grande) ya que puede presentar ambigüedades al compararse con otro kernel. • Aunque normalmente se utilicen kernels del mismo tamaño, para el principio de la señal de compresión basta con utilizar un kernel pequeño, pero al final se debe utilizar uno más grande que abarque todo el desplazamiento de la señal. • No es necesario detectar el desplazamiento de todos los puntos, ya que este no presentará discontinuidades y se puede interpolar. En la siguiente figura se puede observar la señal de deformación calculada, en ella se han representado las zonas válidas en verde y las ambigüedades en azul. En un postprocesado se deberán elegir las zonas válidas y descartar las que presenten ambigüedades. Figura 4.4: Zonas validas y ambigüedades de la señal de desplazamiento acumulado. 37 4.4 Detección de la deformación La principal dificultad es detectar correctamente el desplazamiento de los ecos de la señal de postdeformación respecto a los ecos de la señal de predeformación. Este es uno de los pasos más importantes ya que si se falla aquí no se podrá obtener la deformación del tejido ni representar el elastograma. En el capítulo 3 habíamos hablado de la elasticidad y sus propiedades que venía definida como tensile stress Etensile strain σ ε = = (4.1) Para completar la ecuación que define el módulo de Young definiremos los parámetros de presión y deformación unitaria representados por las letras σ y ε respectivamente. 0 F A E l l σ ε = = ∆ (4.2) De las señales de precompresión y postcompresión, se puede obtener el incremento de tiempo entre reflexiones, el cual nos permiten estimar la deformación. Para simplificar los cálculos y poder medir la deformación estableceremos como constante la velocidad en los diferentes medios, representada por el símbolo c, quedando así la deformación en función del tiempo. ini fin l l l ∆ = − (4.3) l c t = ⋅ (4.4) Combinando las expresiones 4.3 y 4.4 obtenemos la siguiente relación: ( ) fin fin fin ini l c t c t c t t c t ∆ = ⋅ − ⋅ = ⋅ − = ⋅∆ (4.5) Ahora si combinamos las ecuaciones 4.2 y 4.5 obtenemos que el módulo de Young en función de t ∆ . 0 F A E c t l σ ε = = ⋅∆ (4.6) 38 Dado que las señales predeformación y postdeformación han sido digitalizadas deberemos hablar de tiempo discreto (muestras) y no de tiempo continuo (segundos). Conociendo la frecuencia de muestreo podremos pasar fácilmente del dominio discreto al dominio continuo. 1 s s s muestras segundos f T f segundo muestra −     → =         (4.7) S t nT = (4.8) Reagrupando las ecuaciones, podemos establecer una relación entre el dominio discreto y el dominio analógico para obtener t ∆ . s t T n ∆ = ∆ (4.9) Finalmente, substituyendo la ecuación 4.9 en la ecuación 4.6 obtenemos la siguiente ecuación que establece el módulo de Young en función del incremento, o desplazamiento, de muestras. 0 s F A E c T n l σ ε = = ⋅ ⋅∆ (4.10) Para calcular el desplazamiento de muestras partiremos de la señal que hemos obtenido de la interpolación de los puntos válidos, resultado de la función de correlación cruzada. Obtendremos el incremento entre dos puntos sucesivos. En otras palabras, hay que encontrar la pendiente de las rectas de los desplazamientos. Definiremos n ∆ como: ( ) = d deformacion acumulada n dt ∆ (4.11) Como se podemos observar en la ecuación 4.10 para obtener el módulo de Young será necesario controlar varios parámetros como la fuerza aplicada, el área de transductor, la velocidad de propagación… Esto hace que las aplicaciones elastográficas no den el valor del módulo de Young, sino la cantidad de deformación de cada zona. Para este trabajo se han probado dos métodos diferentes para determinar la deformación, uno el uso del factor de correlación y el otro el uso de la función de correlación cruzada. Los resultados son muy parecidos pero la función de correlación cruzada presenta varias mejoras frente al factor de correlación. Definiremos la señal de desplazamiento acumulado como aquella señal que contiene el desplazamiento de un patrón respecto el inicio de la señal. La señal de desplazamiento diferencial es la derivada de la 39 señal de desplazamiento acumulado, es decir, aquella que contiene los incrementos de un punto al anterior. 4.4.1 Factor de correlación Uno de los métodos probados en este trabajo para detectar el desplazamiento temporal ha sido mediante el factor de correlación. Hemos creado un programa que permitía detectar el desplazamiento, pero con un elevado coste computacional y poca tolerancia al ruido. Para grandes rangos de búsqueda ha podido llegar a tardar algunos minutos para una sola A-line. Únicamente se expondrán los algoritmos utilizados pero sin entrar en detalle ni se comentarán los resultados. 4.4.1.1 Fundamentos teóricos del factor de correlación. El coeficiente de correlación de Pearson es un índice estadístico que mide la relación lineal entre dos variables cuantitativas. A diferencia de la covarianza, la correlación de Pearson es independiente de la escala de medida de las variables. El cálculo de correlación lineal se realiza dividiendo la covarianza por el producto de las desviaciones estándar de ambas variables. XY X Y r σ σ σ = ⋅ (4.1) Siendo XY σ la covarianza de (X, Y) y, X σ y Y σ las desviaciones marginales. El valor del índice de correlación varía en el intervalo [-1, +1]: • Si r = 0, no existe relación lineal. • Si r = 1, existe una correlación positiva perfecta. • Si 0 < r < 1, existe una correlación positiva. • Si r = -1, existe una correlación negativa perfecta. • Si -1 < r < 0, existe una correlación negativa. 4.4.1.2 Algoritmo de detección mediante el factor de correlación Basándonos en la teoría del punto anterior hemos desarrollad un algoritmo que permite detectar el desplazamiento por búsqueda de patrones. El algoritmo se detalla a continuación: 1. De la señal precompresión se toma un kernel de duración w k , que tiene como índice de inicio n k y como fin n k + w k . 2. De la señal de postcompresión se toma una muestra de la misma duración que el kernel (w k ) que tenga como inicio m pc , con m pc < n k . 3. Se calcula el factor de correlación entre el kernel y la muestra. 40 4. Se vuelve a tomar una muestra de la señal de postcompresión pero desplazada un instante. 5. Se vuelve a calcular el factor de correlación. 6. Se repiten los pasos 4 y 5 hasta que se llegue al índice de fin de búsqueda establecido. 7. De todos los factores de correlación obtenidos se obtiene la posición del máximo. 8. Se desplaza índice de inicio del kernel de precompresión una muestra. 9. Se repiten los pasos 2 a 7 hasta el final. 10. A partir de todos los máximos del factor de correlación se crea la señal de deformación acumulada. Filtro paso banda centrado en f 0 Leer señal predeformación Filtro paso banda centrado en f 0 Obtener factor de correlación Suma a m pc el desplazamiento de la ventana post Obtener kernel: Duración: w k Inicio: n k Obtener kernel: Duración: w k Inicio: m pc < n k Leer señal postdeformación Suma a n k el desplazamiento de la ventana pre ¿Es m pc mayor que el límite de búsqueda ? Sí No Obtener el máximo y su posición Corregir desplazamientos ¿ Ha recorrido toda la señal? No Sí F igura 4.5 : Algoritmo para detectar el desplazamiento entre señales. 41 Figura 4.6: Kernel y rango de búsqueda para una dimensión. Una vez se ha terminado de correr el kernel de postdeformación sobre la región de búsqueda se obtendrá el máximo de todas las correlaciones, éste indicará en qué punto de la búsqueda se ha producido la máxima coincidencia entre ambas señales. Ya que la ventana ha empezado a buscar unas posiciones antes (m pc < n k ) bastará con restar el desfase inicial entre ambas ventanas para obtener la posición real del máximo que equivale al desplazamiento temporal. Al finalizar todas las correlaciones se obtendrá un vector con todos los desplazamientos temporales. Obsérvese que todos los desplazamientos deberán ir incrementando su valor ya que la señal ha sufrido una compresión, si se redujera el valor del desplazamiento podría tratarse de una ambigüedad. Figura 4.7: Representación de la deformación para el caso ideal sin ambigüedades. 48 49 Capítulo 5. Materiales y métodos. Para probar el funcionamiento de nuestros algoritmos de detección de desplazamiento hemos creado muestras con diferentes estructuras internas. Generalmente estas muestras reciben el nombre de phantom e intentan emular los tejidos humanos y contienen cuerpos esféricos de diferentes diámetros para emular quites o tumores. 5.1 Materiales Los phantoms para ultrasonidos son normalmente de dos tipos. Unos imitan las propiedades acústicas de los tejidos (aproximación a la velocidad del sonido, atenuación media,…). Mientras que los otros intentan simular la apariencia sonográfica del tejido. Estos últimos se suelen utilizar para entrenar la obtención de biopsias o para calibrar los instrumentos de medición. Los phantoms que imitan las propiedades acústicas del sonido han sido fabricados con agar y en el que están suspendidas partículas de grafito, espuma de poliuretano y geles de silicatos de magnesio, y suelen ser utilizados para la evaluación de sonografías o para el estudio de la interacción de los ultrasonidos con los tejidos. Estos phantoms no se suelen utilizar para realizar biopsias ya que suelen ser caros o se suele durar mucho en fabricar. Los phantoms para realizar biopsias suelen ser de fabricación más simple, contienen quistes o masas simuladas y son además ecogénicos o sonoluscentes. Los medios ecogénicos consisten en harina o almidón en una suspensión gelatinosa, polvo de silicium carbide en una suspensión de agar, o simplemente agar. Si están correctamente fabricados, su ecogeneidad simula tejidos parenquimales (excepto agar, que es poco ecogénico y tampoco simula tejido parenquimal); de todas formas, estos pueden ser de producción muy costosa, requieren que se remuevan durante su enfriamiento para que se formen scatterers (dispersores) o que los objetos a biopsia se mantengan suspendidos y requieran materiales que no son de uso muy frecuente (agar, silicium carbide). Medios sonoluscentes o ecotransparentes consisten en gelatinas sin scatterers suspendidos y muy fáciles de preparar. Desafortunadamente son transparentes (a no ser que sean intensamente teñidos), lo que permite que se vean los objetivos o inclusiones. Esto no simula la situación in vivo, y hace que la tarea sea fácilmente realizable para casos de entrenamiento. 50 5.1.1 Concentraciones y formas La elasticidad de la gelatina dependerá de su concentración. Las concentraciones utilizadas en los phantoms B y C las hemos extraído del documento [5] dónde ya se habían caracterizado previamente las gelatinas. Hemos añadido conservante para evitar la descomposición bacteriana y aumentar su vida útil. Tipo de capa Concentración (%) Capa blanda Gelatina 0,05% Nipagin 0,0001% Capa dura Gelatina 0,16% Nipagin 0,0001% Tabla 5.1: Concentraciones de gelatina y conservante. Como recipiente para dar forma a la gelatina hemos elegido un contenedor de plástico impermeable, con la base plana y con tapa hermética para evitar la evaporación del agua. 5.1.2 Elaboración de la muestra En primer lugar hemos estimado la cantidad de agua necesaria para crear una capa. A partir de la cantidad de agua estimada hemos calculado la cantidad de gelatina necesaria. Figura 5.1: Materiales necesarios para la elaboración de phantoms. Para crear la mezcla de gelatina hemos separado un poco del agua y en un recipiente calentamos el resto de agua hasta la ebullición. Disolvimos la gelatina en el agua que previamente habíamos apartado. 51 Finalmente vertimos la gelatina disuelta en el agua caliente y la retiramos de la fuente de calor. Dejamos enfriar la mezcla sin que se endurezca y retiramos la espuma que se haya podido formar. Al verter la segunda una capa de gelatina sobre otra deberemos tomar algunas precauciones para evitar la formación de irregularidades en la superficie. La capa de gelatina sobre la cual vamos a verter la siguiente debe estar perfectamente solidificada y fría. También recomendamos interponer una cuchara u objeto entre el chorro de gelatina líquida y la gelatina solidificada para evitar que se formen agujeros. Como mayor diferencia de temperatura haya entre la gelatina líquida y la solidificada, más limpias o exactas serán las transiciones entre capas. Para la conservación de la gelatina se deberá elegir una zona seca y oscura para evitar la formación de hongos y microorganismos. Se deben eliminar las condensaciones de agua. Si se siguen estas recomendaciones, se puede conservar guardar la gelatina hasta dos semanas. Puede que entre ambas gelatinas haya un intercambio de humedad y así las concentraciones de gelatina-agua variar a lo largo del tiempo. En nuestro caso, hemos observado como el colorante se ha fundido con las diferentes capas; esto nos podría servir como indicador del estado de la gelatina. 5.1.3 Otras consideraciones 5.1.3.1 Conservante para geles. Para aumentar la duración de la gelatina de procedencia animal se propone substituir el agua por agua conservante. Esta agua evita la formación de hongos y bacterias, es una fórmula utilizada en la industria farmacéutica. Agua conservante. Contiene agua purificada con 1% de solución concentrada de p-hidroxibenzoato. [6] Solución concentrada p - hidroxibenzoato Substancia Cantidad P-hidroxibenzoato de metilo (metilparaben – E 206): 8 g 8 g P-hidroxibenoato de propilo (propilparaben – E 216): 1 g 1 g Propilenglicol c.s.p. 100 mL Tabla 5.2: Solución concentrada de p-hidroxibenzoato. Estos ingredientes se pueden encontrar en cualquier farmacia bajo los nombres comerciales de Nipagin (metilparaben) y Nipasol (propiloparaben). Los nombres de las substancias pueden cambiar, en nuestro caso hemos utilizado unicamente metilparaben (Nipagin). 52 5.1.3.2 Dispersores En el capítulo 2 hemos visto que se producirá una reflexión cada vez que haya un cambio de impedancia acústica. En nuestro caso, las gelatinas producirán un eco cada vez que haya un cambio entre capas, de hecho solamente se producen pocas reflexiones. Para aplicaciones más avanzadas a lo mejor es conveniete insertar dispersores (scatters) para que produzcan más ecos y sea más facil identificar determinadas zonas. En función de la longitud de onda de la señal convendrá utilizar unos dispersores u otros. Material Observaciones Semillas de sesamo Producen ecos detectables para frecuencias de 1 MHz. Fibra muciloide hidrofílica No se han producido ecos detectables para frecuencias de 1 MHz. Cambian las propiedades de la gelatina haciendola más blanda y quebradiza. Tabla 5.3: Materiales probados para actuar como dispersores. 5.1.3.3 Inclusiones Producir el phantom por capas es la forma más satisfactoria de producir quistes y masas internos. Una capa de la mezcla se vierte en el contenedor y se deja enfriar hasta que la capa se lo bastante firme como para soportar la inclusión. Si la inclusión flota, se deberá colocar la inclusión en la capa y se verterá el volumen restante de gelatina (sumergida parcialmente en agua caliente para evitar la congelación). Si la inclusión flota, se verterá una capa fina y se introducirá la inclusión, sujetándola. Finalmente se verterá el volumen restante y se dejará enfriar hasta su congelación. Los quistes se simularan con globos llenos de agua, dedos de guantes de exploración, uvas o supositorios de glicerina. Las masas se simularan con trocos de zanahoria, macarrones, aceitunas o trozos de salchichas. Si se refrigera bien, el phantom no sufrirá degeneración bacteriana hasta pasadas unas semanas. 53 5.1.4 Phantoms producidos Figura 5.2: Phantom de dos capas sometido a una pequeña deformación Figura 5.3: Phantom de cinco capas con una sometido a una gran compresión. Como se puede observar en la figura 5.3 la deformación no es igual en todo el phantom, es mayor en la superficie y menor a más profundidad. 54 Figura 5.4: De izq. a dcha.: phantom con dos capas gruesas, phantom con una capa gruesa y phantom con dos capas finas. Figura 5.5: Phantom con dispersores de sésamo y fibras mucilaginosas. 5.2. Montaje experimental 5.2.1. Descripción de los elementos Ordenador con tarjeta NI PCI 5102 – Tarjeta de adquisición de datos La tarjeta de adquisición de datos NI PCI-5102 permite capturar dos canales a 20MS/s a tiempo real con un ancho de banda de 15 MHz. Permite capturar señales con un nivel de ±50 mV hasta ±5 V. 55 Figura 5.6: Ordenador con tarjeta NI PCI 5102. 500PR PULSER-RECEIVER - Emisor de pulsos El 500PR Panametrics es un modelo de alta calidad, emisor-receptor de pulsos ultrasónicos de uso general. Los controles del panel frontal del Modelo 500PR permiten un fácil ajuste de la altura del pulso, aumento de la amortiguación de la onda recibida, y la tasa de repetición de pulsos, así como el modo de funcionamiento (pulso-eco o por transmisión). Un filtro de paso alto seleccionables mediante un interruptor en el panel trasero se puede utilizar para cambiar el tiempo de recuperación del pulso de excitación y de rechazar el ruido de baja frecuencia. Figura 5.7: 500PR vista frontal. 56 Figura 5.8: 500PR vista posterior. V392 – Transductor El transductor receptor V392 es un transductor con un elemento de gran diámetro. Su carcasa exterior tiene un diámetro de 45 mm. Características del transductor V392 Frecuencia 1 MHz Tamaño nominal del elemento 38 mm Point Target Focus Min 6.37 cm; max 19.24 cm Tabla 5.4: Características del transductor V392. Figura 5.9: Transductor V392 57 5.2.2 Montaje y colocación de los elementos Para realizar las medidas se ha realizado el siguiente montaje. 1. El phantom se coloca sobre una superficie plana. 2. Se conecta la salida de la señal del 500PR a la tarjeta de adquisición y la ecosonda a la entrada marcada con T/R. 3. Se conecta la salida de disparo (trigger) del 500PR a la entrada de disparo de la tarjeta de adquisición. 4. Entre ecosonda y phantom pondremos líquido de acoplamiento (glicerina) y moveremos la ecosonda paralelamente a la superficie hasta que las burbujas de aire hayan desaparecido completamente. 5. Ejecutaremos el programa de adquisición de datos en LabVIEW. 5.1 Ajustar parámetros en el programa (cantidad de muestras a capturar y frecuencia de muestreo). Pulsar opción run 6. Visualizar la pantalla la señal y ajustar los parámetros en el 500PR. 5.2 Pulse repetition frequency (PRF): Normalmente al máximo. 5.3 Damping. Normalmente a 5. 5.4 Mode. Emisión y recepción en un mismo punto. Posición inferior 5.5 Pulse height. Nivel del pulso. Hi. 5.6 Gain interruptor: Normalmente a 5. 5.7 Gain potenciómetro: Normalmente a 5. 7. En el programa ajustar los intervalos para guardar la señal ver 5.2.3. 8. Guardar señal predeformación. 9. Con la ecosonda comprimir el phantom y volver a guardar la señal. Figura 5.10: Montaje experimental. 64 Control A. Ruta de acceso a los archivos contenedores de datos, extensión LVM (LabVIEW Measurement File). Se deberá introducir la ruta completa a los archivos guardados con el programa de captura de datos. Figura 5.24: Control de archivos. Control D. Control de retraso o delay con el que se ejecuta el programa para observación de la correlación. Este control de retraso no afecta a los resultados, solamente a la velocidad con la que se ejecuta el programa. Es útil para observar el resultado (sección 2.3) de la función de correlación cruzada entre ambas ventanas (sección 2.1 y 2.2) los resultados de la correlación de las ventanas (sección 2.3). Para un rápido funcionamiento del programa el retraso debe ser 0 (valor por defecto). El indicador Inicio actual indica la posición de inicio de los kernels durante el proceso de cálculo. Figura 5.25: Control del retraso (delay) para visualizar. Control R. Se deberán introducir los índices ente los cuales se van a recortar los resultados. Por defecto ambos valores están a 0, si no se modifican no se mostrarán los resultados. Figura 5.26: Control de recorte Control W. Permite controlar el ancho de los kernels para aplicar la función de correlación cruzada. Muestra el progreso de las correlaciones. Figura 5.27: Control de la ventana de correlación. 65 Control U. Permite controlar los parámetros del detector de picos. Se debe establecer un umbral para que se detecten las zonas de máxima correlación y así poder seleccionar los puntos para la interpolación. Figura 5.28: Control del detector de picos. Ventanas 1.1 y 1.2. Ventanas de visualización para la señal de precompresión y postcompresión. Se muestra la señal sin haber sido filtrada y la señal después de haber sido filtrada por un filtro paso banda centrado en 1 MHz con f L = 900 KHz y f H = 1100 KHz. Figura 5.29: Ventanas de visualización de las señales leídas Los parámetros de los filtros se han establecido por defecto a esas frecuencias. Si se quieren cambiar bastará con acceder al Block Diagram y modificarlos. Los filtros están incluidos en este programa como un VI, por lo tanto, si se cambia la ubicación del archivo también se deberán incluir los filtros. Ventanas 2.1 y 2.2. Ventanas de visualización de los kernels de las señales de precompresión y postcompresión. Figura 5.30: Ventana de visualización de los kernels. 66 Ventanas 2.3 y 2.4. En la ventana 2.3 (derecha de la figura) se muestra el resultado de aplicar la función de correlación cruzada ambos kernels. En la ventana 2.4 se muestra el resultado anterior pero después de haber sido enventanado por una ventana de triangular Hanning para atenuar los picos laterales (lejanos al centro) del resultado de la correlación. Figura 5.31: Ventanas de visualización 2.1 a 2.4. Ventana 3.1. Se visualizan las señales de precompresión y postcompresión superpuestas. Esta ventana es útil para observar la coincidencia de ambas señales y ayuda a seleccionar un ancho de ventana adecuado en el control U. Figura 5.32: Ventana 3.1 Ventana 3.2. El eje horizontal muestra el inicio de los kernels y el eje vertical el valor del máximo del resultado de aplicar la función de correlación cruzada a ambos kernels. Figura 5.33: Ventana 3.2 67 Ventana 3.3. Muestra el desplazamiento respecto del centro de cada uno de los máximos de los resultados de aplicar la función de correlación para cada uno de los kernels. Puede que aparezcan ambigüedades debido a una mala selección del tamaño de la ventana. Se deberán hacer pruebas hasta obtener el menor número de ambigüedades posible. Figura 5.34: Ventana 3.3 Ventana 4.1. Señales de precompresión y postcompresión recortadas con el control de recorte R. Los resultados se actualizarán después de que se vuelva a ejecutar el programa. Las ventanas posteriores solamente mostrarán los resultados para dicho intervalo. Figura 5.35: Ventana 4.1 Ventana 4.2. Valores de los máximos de correlación recortados (ver descripción de la ventana 3.2). Ventana 4.3. Valores de los desplazamientos recortados (ver descripción de la ventana 3.3). Ventana 5.1. Esta ventana muestra lo mismo que las ventanas 4.2 y 4.3 pero de forma superpuesta. Su finalidad es permitir que se seleccione el umbral y el ancho del pico. Figura 5.36: Ventana para determinar el umbral y ancho de pico. 68 Ventanas 5.2 y 5.3. En la primera ventana se pueden ver los puntos seleccionados según los criterios definidos en el control U. En la segunda ventana se mostrará la señal interpolada. Uno de los requisitos para que la elastografía se obtenga correctamente es que el valor de un punto no debe ser mayor que el anterior. En caso de que se observaran puntos que no cumplen esta condición se deberán modificar los valores del control de umbral U. Figura 5.37: Puntos de interpolación incorrectos debido a una mala selección del umbral y el ancho de los picos. Como se puede observar en la figura 41 los puntos de interpolación no cumplen la condición que hemos comentado anteriormente. Los valores empleados para seleccionar esos puntos han sido 0 para el umbral y 3 muestras para el ancho de los picos. Mediante el cursor de la ventana 5.1 obtenemos el valor de Y del menor pico que es de 13,62, y el ancho del pico son unas 300 muestras. Ahora estableceremos como umbral unos valores ligeramente inferiores para que se detecten correctamente los desplazamientos y se descarten las ambigüedades o espúrreos. Unos buenos valores serían un umbral de 12 y un ancho del pico de 200 muestras. Figura 5.38: Resultado de seleccionar umbrales y anchos de picos adecuados. Figura 5.39: Ventana 5.3. Señal interpolada Todos los segmentos deberán tener el mismo signo de la pendiente. El signo de la pendiente puede cambiar dependiendo si se aplica la función de correlación de un conjunto A sobre B o B sobre A. 69 Ventana 5.4. Muestra el resultado de derivar la señal mostrada en la ventana 5.3. Los valores cercanos a 0 representan poca deformación mientras que valores mayores en módulo representan deformaciones mayores. Figura 5.40: Derivada de los desplazamientos acumulados. Ventana 6. Última ventana del programa, elastograma o mapa de deformación de las señales cargadas al principio. Esta ventana representa los valores de obtenidos en la ventana 5.4 según una escala de color. Para que se puedan observar bien las diferencias se deberán ajustar los marcadores de la escala de color. Figura 5.41: Elastograma unidimensional. 70 5.5 Simulación numérica Para comprobar que los resultados obtenidos a partir de señales reales son válidos y si influye la presencia de ruido en los resultados. Se ha creado un modelo virtual y se ha estudiado la propagación de la señal en este utilizando la técnica de las diferencias finitas en el dominio del tiempo (FDTD). El programa de simulación está realizado en MATLAB y consiste en un script llamado simulacion.m y una función llamada fdtd2d.m. 5.5.1 FDTD (Finite Difference Time Domain) Las diferencias finitas en el dominio del tiempo (FDTD) es una técnica de modelización computacional muy utilizada. Se considera fácil de entender y fácil de implementar en software. Dado que es un método que trabaja en el dominio temporal, las soluciones pueden abarcar un amplio rango de frecuencias con una sola ejecución de simulación. El método FDTD pertenece a la clase general del modelado numérico basado en diferencias finitas de mallas. Para ello se discretizan las ecuaciones de onda y se traducen a ecuaciones de diferencias finitas. Posteriormente se resuelven todas las ecuaciones implicadas, repitiendo y actualizando, para cada instante hasta que se consiguen los resultados deseados 5.5.2 Programa de simulación fdtd2.m simulacion.m Calcula y actualiza mallas de velocidad y presión Envía parámetros Calcula diferenciales Devuelve señales y fotogramas si se piden. Enventana señales y calcula FCC para cada ventana Desplazamiento acumlado Seleccionamos puntos Desplazamiento diferencial ELASTOGRAMA Figura 5.42 : Funcionamiento simplificado del programa de simulación. 71 El script simulacion.m está formado por tres partes. Las dos primeras (parámetros de simulación y parámetros de procesado) serán las que se deban modificar para controlar las señales y su postprocesado. Los parámetros y cómo influyen en las señales se comentan en el script simulacion.m y en la función fdtd2d.m Este programa obtiene dos señales, una de precompresión y otra de postcompresión. La señal de precompresión se obtiene a partir de la simulación de una estructura bidimensional de tres capas donde se propaga un pulso ultrasónico, produciéndose una reflexión con cada cambio de impedancia. La obtención de la señal de postcompresión sigue el mismo principio que la señal de precompresión pero sufre una compresión, que se indica en los parámetros que se envían a la función fdtd2d.m. Una vez que se han obtenido las señales se procese al postprocesado de la misma forma que lo realiza el programa en LabVIEW, utilizando kernels y aplicando la función de correlación cruzada. Se obtiene la señal de desplazamientos acumulados y de esta se interpolan los puntos (previa selección) para realizar la primera derivada y obtener la señal de desplazamiento diferencial. Finalmente se representa la señal de desplazamiento diferencial según un mapa de colores. Figura 5.43: Ventana de presentación de los resultados del programa de postprocesado en MATLAB. 72 73 Capítulo 6. Resultados y discusión. 6.1 Creación de phantoms. Durante la elaboración del proyecto se han creado varios phantoms siguiendo los pasos descritos en el punto 5.1.2. Visualmente se han podido diferenciar las diferentes capas por su coloración, las capas blandas de un color más claro y las capas duras de un color más oscuro. Uno de los phantoms ha tenido que ser desechado ya que durante el proceso de elaboración la gelatina aún caliente ha derretido la capa inferior creando un agujero en esta. También se ha han podido diferenciar las diferentes elasticidades de las capas mediante el tacto. En los primeros phantoms creados se añadió colorante alimentario a las capas duras para poder diferenciar mejor las diferentes capas. Al cabo de una semana se ha observado como el colorante ha coloreado las capas blandas sin colorante. La coloración de las capas podría indicar un intercambio de agua entre las diferentes capas lo que podría tener como consecuencia un cambio de la elasticidad. A pesar de ello se continuaban detectando las diferentes capas mediante el método elastográfico. Se ha creado un phantom con dispersores a base de fibras mucilaginosas. La concentración de gelatina ha sido la misma que para los demás phantoms, las cantidades de fibras incluidas en el phantom han sido las recomendadas en la sección 5.1.1. En este phantom se ha observado en mediante el tacto que la diferencia de elasticidad no era tan grande como en las gelatinas sin fibras mucilaginosas. La superficie de contacto del phantom con el transductor se quebró durante las pruebas. Debido a los colorantes que se adjuntaban como excipientes a las fibras mucilaginosas todas las capas del phantom adquirieron un color marrón oscuro haciendo que fuera más difícil diferenciar las capas visualmente (ver figura 6.2). Después de hacer varias medidas y observar que no se obtenían más ecos que sin fibras, se decidió añadir dispersores más grandes (semillas de sésamo). Se volvió a fundir la gelatina con las fibras mucilaginosas y se añadieron las semillas de sésamo. En este segundo phantom se observaron más ecos debido a las reflexiones de las semillas pero al final se desestimó la utilización de las semillas ya que no se creyó que no era adecuado para exploraciones en una dimensión. Para futuros trabajos donde se realicen exploraciones en más de una dimensión y a mayor frecuencia sí que se pueden añadir dispersores. En cuanto a la vida útil de los phantoms, los primeros se realizaron sin conservantes y su duración fue de una semana en nevera. Pasado este tiempo se observaron síntomas de descomposición por la acción bacteriana. Los últimos phantoms creados, con conservantes, han llegado a durar hasta tres semanas fuera de la nevera sin presentar signos de descomposición. En algunos el agua condensada disolvió la gelatina pero tras un lavado se pudo volver a fundir la gelatina y crear de nuevo el phantom. Finalmente, podemos concluir que la utilización de gelatinas para la creación de phantoms es una buena opción ya que permiten la creación de phantoms según las necesidades, de manera rápida, sencilla y económica. Los diseños también pueden ser variados, desde gelatinas sencillas en capas a gelatinas con inclusiones esféricas o de cualquier otra forma. El precio del phantom ha sido de aproximadamente de 6 80 Figura 6.9: Resultados parciales y elastograma para la simulación 2. Simulación 2 – Deformación 1 % 81 Figura 6.10: Resultados parciales y elastograma para la simulación 3. Simulación 3 – Deformación 2 % 82 Figura 6.11: Resultados parciales y elastograma para la simulación 4. Simulación 4 – Deformación 4 % - Simulación de control 83 6.4 Postprocesado de datos reales Para comprobar que el programa de postprocesado detecta correctamente los desplazamientos realizaremos la comparación entre los resultados del programa y los resultados mediante el recuento manual. A continuación también se estudiará la influencia del tamaño de la ventana de correlación y se justificará el tamaño elegido. Finalmente se expondrán los elastogramas obtenidos. 6.4.1 Comparación de los resultados del programa con resultados de recuento manual Para demostrar el correcto funcionamiento del programa de postprocesado de datos compararemos los resultados de la obtención del retraso acumulado obtenidos mediante el programa con los resultados obtenidos mediante recuento manual. Una vez cargadas las señales de predeformación y postdeformación de la medida 1.1.1 ejecutaremos el programa, dejaremos la ventana que viene por defecto con unas 600 muestras de ancho. Figura 6.12: Señales de precompresión y postcompresión superpuestas. En la señal de la figura anterior se observan cuatro zonas, la primera zona corresponde al pulso emitido y las otras tres a los ecos. Como es de esperar los pulsos emitidos se solapan perfectamente (ver principio de la señal) y la señal de postcompresión está desplazada respecto a la de precompresión. La señal obtenida a partir del valor máximo de todas las correlaciones también presenta cuatro picos que corresponde a las zonas descritas anteriormente. Figura 6.13: Valores de las correlaciones. 84 La siguiente señal que se va a analizar es la de los desplazamientos acumulados, resultado del desplazamiento de los picos de la función de correlación cruzada respecto del centro. Figura 6.14: Desplazamiento acumulado. En la señal de desplazamiento acumulado se observan varias zonas planas localizadas en las regiones que hemos comentado para las señales de predeformación y postdeformación y el valor de correlación. Si solapamos ambas señales se observa que algunas zonas planas corresponden con los picos de correlación. Por otra parte en la señal de desplazamientos acumulados presenta varios picos positivos y negativos, esto es debido a que los patrones de la señal se repiten y puede que tengan una correlación mayor para zonas lejanas al desplazamiento real. De antemano también se pueden descartar todos los picos positivos ya que si la señal se ha comprimido los picos necesariamente deben ser negativos, un pico positivo indica expansión en la señal o retraso. En el caso de la simulación la señal de correlación era de signo positivo, esto es debido al proceso de obtención de la posición del máximo (ver resultados de las simulaciones). Figura 6.15: Señales solapadas. Como se había comentado en su momento, una de las propiedades de la FCC como resultado cero si se aplica a dos conjuntos de valores aleatorios (ruido blanco). Estos valores cercanos a cero nos indicarán la validez de los retrasos calculados. Por tanto, descartaremos los valores de los retrasos que estén en zonas con un valor máximo de correlación igual o cercano a cero y tomaremos como válido aquellos valores que estén por encima de un umbral de valor del máximo de correlación. En esta sección tomaremos como válidos los valores de retraso en aquellas zonas donde están los ecos de las señales y los picos de la señal de correlación. Mediante cursores obtendremos los valores de los 85 retrasos en las zonas donde existen máxima correlación y posteriormente los comparemos con la medida manual. Figura 6.16: Señales solapadas con cursores. Retrasos calculados por el programa Índice Valor Cursor A 105 -0.5 Cursor B 3208 -214.5 Cursor C 4906 -217.5 Cursor D 7484 -244.5 Tabla 6.3: Retrasos calculados por el programa. Para saber si los valores son correctos compararemos los valores obtenidos de forma programática (señal retrasos) con la diferencia de índices de patrones reconocibles de las señales precompresión y postcompresión. Figura 6.17: Señal completa con los cursores Al igual que de forma programática hemos establecido cuatro zonas donde se pueden identificar los ecos. Procedemos a calcular la diferencia entre el índice de llegada del eco de precompresión (cursor en mayúscula) y el de postcompresión (cursor en minúscula). 86 Figura 6.18: Zoom cursores A y a Figura 6.19: Zoom cursores B y b Figura 6.20: Zoom cursores C y c Figura 6.21: Zoom cursores D y d 87 Retrasos obtenidos de forma manual Cursor A Cursor a Cursor B Cursor b Cursor C Cursor c Cursor D Cursor d Índice 278 278 3534 3322 5293 5080 7866 7623 Diferencia a - B = 0 b - B = -212 c - C = -213 d – D = -243 Tabla 6.4: Retrasos calculados de forma manual. Comparación de los resultados Programática Manual Diferencia Zona A -0.5 0 0.5 Zona B -214.5 -212 2.5 Zona C -217.5 -213 4.5 Zona D -244.5 -243 1.5 Tabla 6.5: Comparación de los retrasos obtenidos por el programa con los retrasos calculados de forma manual. Las diferencias entre los resultados obtenidos de forma programática y de forma manual es muy pequeña, 4.5 muestras a una frecuencia de muestreo de 100 MHz equivale a 40 ns. Al ser las diferencias tan pequeñas podemos decir que los resultados son válidos y los errores producidos son debidos a la colocación de los cursores. 6.4.2 Cálculo de la velocidad de propagación. Un aspecto muy interesante es obtener la velocidad de propagación de las ondas acústicas ya que la finalidad del phantom es que sea lo más parecido posible a la de los tejidos humanos. Para ello hemos decido hacer las medidas en un phantom sencillo, en este caso la gelatina con la capa gruesa, ya que podremos identificar mejor las reflexiones. En un primer lugar calculamos el tiempo transcurrido en muestras desde el envío del pulso hasta el fondo del phantom. Figura 6.22: Identificación del inicio y el final de la señal. 88 Valores de los cursores Inicio Final Índice 5052 12789 Diferencia Final – Inicio = 7737 Tabla 6.6: Duración del recorrido acústico en muestras. Considerando el que los ecos han viajado dos veces la distancia hasta el obstáculo, la frecuencia de muestreo y las dimensiones de la gelatina obtenemos la siguiente relación: 2 s d c f muestras ⋅ =⋅ (6.1) Dado el espesor del phantom que es de 6 cm y la cantidad de muestras podemos calcular la velocidad del sonido que es de 1550 m·s -1 . Comparando con la tabla 1 podemos decir que la velocidad de propagación es muy cercana a la velocidad de propagación en tejidos humanos. 6.4.3 Elastogramas obtenidos a partir de datos reales Casi todos los elastogramas obtenidos a partir de señales reales se asemejan al phantom del cual se han obtenido las señales. En aquellos casos que ha sido necesaria una ganancia mayor de señal los elastogramas presentaban algunas incoherencias, como zonas más duras. Este fenómeno ya se ha comentado en el apartado de ruido eléctrico introducido por el 500PR. El desplazamiento de que se observa entre el elastograma calculado y la imagen real (fotografía) es debido al recorte de la señal o al escalado de la fotografía. En algunos casos se observan franjas de color del mismo tono, p. ej. medida 1.2.2.. Esto es debido a que el programa selecciona de forma automática los puntos de interpolación superiores a un umbral en los picos de los máximos de correlación. Puede que existan dos picos de máximos de correlación muy cercanos, dado a la presencia de señal, y que por ello se seleccionen dos puntos cercanos. Estos puntos son válidos ya que presentan las mismas características de elasticidad que los puntos cercanos. A mayor cantidad de puntos mayor cantidad de información. Por último comentar que en la zona inferior de los elastogramas se detectan zonas donde la deformación es menor, p. ej. medida 1.1.2, esto se debe a la presencia de un campo de deformaciones, donde la deformación es mayor en la superficie (zona cercana al transductor compresor) que en las regiones más profundas. En la siguiente tabla se recogen a modo de resumen los parámetros más importantes para el postprocesado. 89 Resumen de deformaciones, tamaños de kernel y ganancias. Medida Deformación Tamaño kernel Ganancia 1.1.1 2 mm 600 muestras 4 1.1.2 4 mm 600 muestras 4 1.2.1 2 mm 600 muestras 4 1.2.2 4 mm 600 muestras 4 2.1.1 2 mm 600 muestras 4 2.1.2 4 mm 600 muestras 4 2.2.1 2 mm 600 muestras 4 2.2.2 4 mm 800 muestras 4 3.1.1 2 mm 800 muestras 7 3.1.2 4 mm 800 muestras 7 3.2.1 2 mm 400 muestras 7 3.2.2 4 mm 700 muestras 7 4.1 4 mm 600 muestras 7 4.2 7 mm 1000 muestras 7 4.3 12 mm 1000 muestras 7 Tabla 6.7: Resumen de deformaciones, tamaños de kernel y ganancias. Para cada phantom se han capturados las señales de predeformación (deformación 0 mm) y deformación de 2 y 4 mm. En el último caso, medida 4, se acoplaron dos phantoms para formar uno más grande. En ese caso fue necesaria una compresión mayor para que la deformación llegara al fondo del phantom. Una compresión excesiva hace que se pierda correlación de las señales para las capas superiores donde la deformación es mayor. También hace necesario el uso de ventanas mayores con lo que se pierde resolución espacial (1000 muestras). En las medias 3.1.1, 3.1.2 y 3.2.1 se puede observar que el principio del elastograma (parte superior) muestra valores de poca deformación. Si nos fijamos, el parámetros que tienen las tres medidas en común es que la ganancia es de 7 en lugar de 4 como en las otras medidas. Al ser la ganancia más alta, también se introduce más ruido eléctrico, dando como resultado una información incorrecta. Uno de los aspectos a resaltar del postprocesado es la detección automática de puntos para la interpolación, ya que solamente basta con fijar un umbral en el nivel de correlación. 96 Medida 2.2.1: Gelatina capas gruesas – Base ancha arriba y compresión 2 mm Tipo de señal Nombre del archivo Señal predeformación: A1 - fs 100 M - dur 50 k - def 0mm - gain 4 - damping 0_001.lvm Señal postdeformación: B1 - fs 100 M - dur 50 k - def 2 mm - gain 4 - damping 0_001.lvm Parámetros y valores del postprocesado Parámetro Valor Intervalos de recorte (índices) Inicio Fin 5000 15000 Compresión aplicada 2 mm Ancho de la ventana de correlación 600 Umbral de detección de picos 14 Ancho de los picos 200 muestras Figura 6.29: Comparación elastograma e imagen real para la medida 2.2.1. 97 Medida 2.2.2: Gelatina capas gruesas – Base ancha arriba y compresión 4 mm Tipo de señal Nombre del archivo Señal predeformación: A1 - fs 100 M - dur 50 k - def 0mm - gain 4 - damping 0_001.lvm Señal postdeformación: B1 - fs 100 M - dur 50 k - def 2 mm - gain 4 - damping 0_001.lvm Parámetros y valores del postprocesado Parámetro Valor Intervalos de recorte (índices) Inicio Fin 5000 15000 Compresión aplicada 4 mm Ancho de la ventana de correlación 800 Umbral de detección de picos 12 Ancho de los picos 50 muestras Figura 6.30: Comparación elastograma e imagen real para la medida 2.2.2. 98 Medida 3.1.1: Gelatina capas finas – Base ancha abajo y compresión 2 mm Tipo de señal Nombre del archivo Señal predeformación: A1 - fs 100 M - dur 50 k - def 0 mm - gain 7 - damping 0_001.lvm Señal postdeformación: B1 - fs 100 M - dur 50 k - def 2mm - gain 7 - damping 0_001.lvm Parámetros y valores del postprocesado Parámetro Valor Intervalos de recorte (índices) Inicio Fin 5000 14000 Compresión aplicada 2 mm Ancho de la ventana de correlación 800 Umbral de detección de picos 4 Ancho de los picos 200 muestras Figura 6.31: Comparación elastograma e imagen real para la medida 3.1.1. 99 Medida 3.1.2: Gelatina capas finas – Base ancha abajo y compresión 4 mm Tipo de señal Nombre del archivo Señal predeformación: A1 - fs 100 M - dur 50 k - def 0 mm - gain 7 - damping 0_001.lvm Señal postdeformación: C1 - fs 100 M - dur 50 k - def 4 mm - gain 7 - damping 0_001.lvm Parámetros y valores del postprocesado Parámetro Valor Intervalos de recorte (índices) Inicio Fin 5000 14000 Compresión aplicada 4 mm Ancho de la ventana de correlación 800 Umbral de detección de picos 4 Ancho de los picos 200 muestras Figura 6.32: Comparación elastograma e imagen real para la medida 3.1.2. 100 Medida 3.2.1: Gelatina capas finas – Base ancha arriba y compresión 2 mm Tipo de señal Nombre del archivo Señal predeformación: A1 - fs 100 M - dur 50 k - def 0 mm - gain 7 - damping 0_001.lvm Señal postdeformación: B1 - fs 100 M - dur 50 k - def 2 mm - gain 7 - damping 0_001.lvm Parámetros y valores del postprocesado Parámetro Valor Intervalos de recorte (índices) Inicio Fin 5000 14000 Compresión aplicada 2 mm Ancho de la ventana de correlación 400 Umbral de detección de picos 4 Ancho de los picos 200 muestras Figura 6.33: Comparación elastograma e imagen real para la medida 3.2.1. 101 Medida 3.2.2: Gelatina capas finas – Base ancha abajo y compresión 4 mm Tipo de señal Nombre del archivo Señal predeformación: A1 - fs 100 M - dur 50 k - def 0 mm - gain 7 - damping 0_001.lvm Señal postdeformación: C1 - fs 100 M - dur 50 k - def 4 mm - gain 7 - damping 0_001.lvm Parámetros y valores del postprocesado Parámetro Valor Intervalos de recorte (índices) Inicio Fin 5000 14000 Compresión aplicada 4 mm Ancho de la ventana de correlación 700 Umbral de detección de picos 100 Ancho de los picos 200 muestras Figura 6.34: Comparación elastograma e imagen real para la medida 3.2.2. 102 Medida 4.1: Doble gelatina – Compresión 4 mm Tipo de señal Nombre del archivo Señal predeformación: A1 - fs 100 M - dur 50 k - def 0 mm - gain 7 - damping 0_001.lvm Señal postdeformación: B1 - fs 100 M - dur 50 k - def 4 mm - gain 7 - damping 0_001.lvm Parámetros y valores del postprocesado Parámetro Valor Intervalos de recorte (índices) Inicio Fin 5000 21000 Compresión aplicada 4 mm Ancho de la ventana de correlación 1000 Umbral de detección de picos 50 Ancho de los picos 200 muestras Figura 6.35: Comparación elastograma e imagen real para la medida 4.1. 103 Medida 4.2: Doble gelatina – C ompresión 7 mm Tipo de señal Nombre del archivo Señal predeformación: A1 - fs 100 M - dur 50 k - def 0 mm - gain 7 - damping 0_001.lvm Señal postdeformación: C1 - fs 100 M - dur 50 k - def 7 mm - gain 7 - damping 0_001.lvm Parámetros y valores del postprocesado Parámetro Valor Intervalos de recorte (índices) Inicio Fin 5000 21000 Compresión aplicada 7 mm Ancho de la ventana de correlación 1000 Umbral de detección de picos 100 Ancho de los picos 300 muestras Figura 6.36: Comparación elastograma e imagen real para la medida 4.2. 104 Medida 4.3: Doble gelatina – Compresión 12 mm Tipo de señal Nombre del archivo Señal predeformación: A1 - fs 100 M - dur 50 k - def 0 mm - gain 7 - damping 0_001.lvm Señal postdeformación: D1 - fs 100 M - dur 50 k - def 12 mm - gain 7 - damping 0_001.lvm Parámetros y valores del postprocesado Parámetro Valor Intervalos de recorte (índices) Inicio Fin 5000 21000 Compresión aplicada 12 mm Ancho de la ventana de correlación 1000 Umbral de detección de picos 40 Ancho de los picos 200 muestras Figura 6.37: Comparación elastograma e imagen real para la medida 4.3. 105 6.4.4 Influencia del tamaño del kernel en la determinación de los retrasos En esta sección se detallarán de forma breve la influencia del tamaño del kernel en la determinación de los retrasos para la obtención de la imagen elastográfica. Resumen de los tamaños de kernel Deformación 2mm Deformación 4 mm 100 100 200 200 400 400 600 600 800 800 1000 1000 1600 1600 3200 3200 6400 6400 1000 1000 Tabla 6.8: Resumen de los tamaños de kernel. En los casos que se utilizan kernels pequeños en la señal de desplazamiento acumulado aparecen muchos picos y ambigüedades. A medida que el tamaño del kernel aumenta van apareciendo zonas planas o mesetas donde el desplazamiento se mantiene constante. A partir de tamaños de unas 550 muestras, para este trabajo, se puede decir que el tamaño de las mesetas en la señal de desplazamiento acumulado es suficiente para poder seleccionar de forma automática, a partir del nivel de correlación, los puntos para la interpolación. Para tamaños muy grandes de kernel se obtienen mesetas muy anchas con el mismo nivel de desplazamiento. Esto es debido a que las ventanas de correlación abarcan un mismo patrón, con mucho nivel, y si las dos van desplazándose el desplazamiento permanece constante. Esta meseta se extenderá hasta que se abarque con el kernel otro patrón. Finalmente debemos decir que la señal de desplazamiento acumulado debe ser considerada conjuntamente con la señal de máximos de correlación. La señal de máximos de correlación nos indicará la validez de las zonas de desplazamiento acumulado. El tamaño de la ventana debe ser suficiente para que abarquen dos patrones reconocibles, con la suficiente información para que se puedan considerar únicos y se puedan correlar ambas ventas. 112 [13] E. KONOFAGOU ELISA; OPHIR JONATHAN; A. KROUSKOP THOMAS, S. GARRA BRIAN. Elastography: from theory to clinical applications. 2003 Summer Bioengineering Conference, June 25-29, Sonesta Beach Resort in Key Biscayne, Florida [14] DANIEL T. GINAT; HUNG GILBERTO; R. GARDNER THOMAS; E. KONOFAGOU ELISA. High resolution ultrasound elastography of articular cartilage in vitro. Department of Biomedical Engineering, Columbia University, New York. [15] ABBAS SAMANI; PLEWES DONALD. An inverse problem solution for measuring the elastic modulus of intact ex vivo breast tissue tumours. University of Western Ontario, London, Ontario, N6A 5C1, Canada. April 2006 [16] KYUNG-CHAN KIM; SE BAEK OH; SOO HYUN KIM; YOON KEUN KWAK. Design of a signal processing algorithm for error-minimized optical triangulation displacement sensors. Department of Mechanical Engineering, Korea Advanced Institute of Science and Technology. [17] CHUXIONG DING; JING BAI. Reduction of the decorrelation effect due to tissue lateral displacement by 2-D spatial comprehensive correlation in elastography. Institute of Biomedical Engineering, Department of Electrical Engineering, Tsinghua University [18] YOUSUKE SATO; MICHITSUGU MORI; YASUSHI TAKEDA; KOICHI HISHIDA; MASANOBU MAEDA. Signal processing for advanced correlation ultrasonic velocity profiler. Third International Symposium on Ultrasonic Doppler Methods for Fluid Mechanics and Fluid Engineering EPFL, Lausanne, Switzerland, September 9 - 11, 2002 [19] SOSA-CABRERAA DARIO; TRISTAN-VEGAB ANTONIO; VEGAS-SANCHEZ-FERREROB GONZALO. GONZALEZ-FERNANDEZA JAVIER; GOMEZ-DENIZA LUIS; ALBERLA LOPEZB CARLOS; RUIZ-ALZOLACA JUAN. A new approach to elastography using a modified demons. Center for Technology in Medicine, Dept. Señales y Comunicaciones. [20] CARL D. LATINO; NIVEL LERVAD ANDERSEN; FRANDS VOSS. Detection of time arrival of ultrasonic pulses. Journal of Physics: Conference Series 52 (2006) 14–26 [21] J. GONZÁLEZ, J.M. CANINO, L. GÓMEZ. Una experiencia docente para instrumentación electrónica aplicada a la técnica médica de elastografía por ultrasonido. CTM (Centro de Tecnología Médica) de Universidad de Las Palmas de G.C (ULPGC). España. [22] M. M. ARTIGAO; E. RUBIO; S. DÍEZ; V. CASELLES. Técnicas de correlación cruzada en el registrado de imágenes de diferente resolución espacial. XI Congreso Nacional de Teledetección, 21-23 septiembre 2005. Puerto de la Cruz. Tenerife. 113 [23] MARTÍNEZ PATRICIA; AZUAGA MARCELO. Medición del módulo de elasticidad de Young. Laboratorio IV - Dpto. de Física - UBA 1997 [24] MICHAEL F. INSANA. Elasticity imaging. University of Illinois at Urbana-ChampaignUrbana, Illinois [25] JEREMY C HEBDEN; BEN D PRICE; ADAM P GIBSON; GARY ROYLE. A soft deformable tissueequivalent phantom for diffuse optical tomography. Department of Medical Physics & Bioengineering, University College London. Agosto 2006. [26] TOSHIO KONDO; MITUYOSHI KITATUJI. New Tissue Mimicking Materials for Ultrasound. Medical Div. Takiron [27] BEVERLY STERN. The Basic Concepts of Diagnostic Ultrasound. [28] JAMAL RAHMAN. LabVIEW: applications and solutions. Prentice Hall. ISBN: 0130964239 [29] CHUGANI MAHESH L.. LabVIEW signal processing. Prentice Hall PTR. ISBN: 0139724494 [30] ALMENAR TERRÉ, VICENÇ. Apuntes de sistemas lineales. Editorial Universidad Politécnica de Valencia. 2000. ISBN: 8477219087 [31] KOCI STEFAN; CHAPMAN & HALL. Ultrasonic measurements and technologie. Chapman & Hall. ISBN: 0412638509 [32] BERTRÁN ALBERTÍ, EDUARD. Procesado digital de señales. Fundamentos para comunicaciones y control II. Edicions de la Universitat Politècnica de Catalunya. Año 2006. PÁGINAS WEB [32] Emisores-receptores controlados por computadora 5800PR, 5900PR, 5910PR <http://www.olympus-ims.com/es/5800pr/> [33] CIRS Incorporated <http://www.cirsinc.com/> [34] DR. JOSEP WOO. History of ultrasound obstetrics and gynecology, Part 1. <http://www.ob-ultrasound.net/history.html> [35] Wikipedia EN: Stress (mechanics) http://en.wikipedia.org/wiki/Stress_%28mechanics%29 114 115 ANEXO I: simulacion.m %Simulación de deformación en gelatinas % %Autor: Javier Redondo Pastor %Coautor: Daniel Serra Freese %Abril 2010. PFC Obtención de imágenes elastográficas en 1D. % %Para más información consulta el PFC y consulta el fichero. close all %Cerramos todos los plots clc %Limpiamos linea de comandos clear all %Borramos variables %PARÁMETROS SIMULACIÓN %Generamos las señales de pre y postdeformación [h,w,dt,dh,imagenpre]=fdtd2d(100,80,0,0,'PREDEF'); %Genera la señal de precompresión [h2,w,dt,dh2,imagenpost]=fdtd2d(100,80,2,0,'POSTDEF'); %Genera la señal de postcompresión %OUTPUTS: %[h,w,dt,dh,imagen]=fdtd2d(mtt,maxttt,def,captureMovie,nombre) %h -> alto del material %w -> ancho %dh -> differencial espacial %imagen -> structure con frames %INPUTS: %maxtt -> control bucle interno %maxttt -> control bucle externo, número de frames %def -> deformación en tanto porciento p. ej. 5 (5%) %captureMovie -> Boolean 0 = No se captura; 1 = Se captura %nombre -> Nombre para identificar graficos %PARÁMETROS POSTPROCESADO step = 10; %Fijamos pasos entre correlaciones winSize = 600; %Tamaño del kernel %Guardamos datos save datosguardados %Lo guarda todo en archivo datosguardados.mat load datosguardados %Carga datosguardados.mat %DE AQUÍ EN ADELANTE NO SE DEBE MODIFICAR %Postprocesado t=(1:1000)*dt; %Vector de tiempos x=exp(-1d12*(t-.2d-5).^2).*sin(2*pi*t*2e6); %Señal excitadora %Convolucionamos las señales de pre y post con la señal excitadora captura_pre=conv(x,h); %Convoluciona captura_post=conv(x,h2); %Convoluciona %Obtenemos una señal más ancha y con más información. El resultado %debe ser muy similar si la señal excitadora es del mismo %tipo que el pulso mandando con fdtd2d. 116 for ii=1:step:length(captura_pre)-1-winSize; %Enventana, correla y plot x=captura_pre(ii:ii+winSize); %Enventana pre y=captura_post(ii:ii+winSize); %Enventana post currentCorr=xcorr(x,y); %Antes calculin, correlación actual plot((-winSize:winSize),currentCorr); title(['Correlograma (Inicio actual: ' num2str(ii) ' de ' num2str(length(captura_pre)+winSize) ')']) %Ponemos título pause(.001) %dur{currentCorr} = dur{x} + dur{y} - 1 [corrAmount(ii), desplazamiento(ii)]=max(currentCorr); %Revisar si basura(ii) es la cantidad de correlación %desplazamiento -> vector con los desplazamientos respecto del centro end %Representamos desplazamiento y máximo de correlación subplot(2,2,1) plot((1:step:length(desplazamiento)),desplazamiento(1:step:length(desplazamient o))-winSize); title(['Desplazamiento acumulado']) %Ponemos título subplot(2,2,2) plot((1:step:length(desplazamiento)),10*log(corrAmount(1:step:length(desplazami ento)))) title(['Máximo de correlación dB']) %Ponemos título %dB para acercar valores lejanos disp('PULSA UNA TECLA Y SELECCIONA EL CENTRO DE LAS MESETAS.') disp('El elastograma variará según la zona que pulses, las pendientes') disp('calculadas no serán las mismas.') pause %Se espera hasta que se pulse una tecla [X, Y]=ginput(3); %Pedimos que introduzca 3 puntos con el ratón Xprim(1)=0; Yprim(1)=0; %Para que haya un punto en (0,0) for i=2:1:4 %Insertamos todos los puntos en un único vector [0 X] y [0 Y] Xprim(i)=X(i-1); Yprim(i)=Y(i-1); end interpol=interp1(Xprim,Yprim,1:1:max(Xprim)); %Interpolamos para el siguiente paso subplot(2,2,3) plot(interpol) %Representa todos los puntos interpolados hold on plot(Xprim, Yprim,'go') %Representa los puntos originales title(['Puntos interpolados']) %Ponemos título hold off subplot(2,2,4) diffinterpol=diff(interpol); %Cálculo de la diferencia (pendiente segmento) plot(diffinterpol) title(['Desplazamiento diferencial']) %Ponemos título figure imagesc (diffinterpol'); colorbar; %Representa el elastograma para compararlo. title(['Elastograma']) %Ponemos título disp('') 117 disp('Si quieres ver la pelicula escribe: movie(elnombre)') % hold on; plot(diff(Y(1:2))/diff(X(1:2))*(1:10000)+Y(1)- X(1)*diff(Y(1:2))/diff(X(1:2)),'r') % % X=X([1 3]);Y=Y([1 3]); % hold on; plot(diff(Y(1:2))/diff(X(1:2))*(1:10000)+Y(1)- X(1)*diff(Y(1:2))/diff(X(1:2)),'r') % % [X,Y]=ginput(3); % % hold on; plot(diff(Y(1:2))/diff(X(1:2))*(1:10000)+Y(1)- X(1)*diff(Y(1:2))/diff(X(1:2)),'r') % % X=X([1 3]);Y=Y([1 3]); % hold on; plot(diff(Y(1:2))/diff(X(1:2))*(1:10000)+Y(1)- X(1)*diff(Y(1:2))/diff(X(1:2)),'r') 118 119 ANEXO 2: fdtd2d.m function [h,w,dt,dh,framesCollector]=fdtd2d(maxtt,maxttt,niveldecompresion,captureMovie, name) %Define dimensiones de capas de materiales capa1=.01; %Metros capa2=.02; capa3=.025; j = 1; %Indice para los frames dh=(.25d-3)/4; %Resolución espacial 0.0625 milimetros nx=round(.0014/dh); %Indices x, a lo horizontal ANCHOOOO ny1=round(capa1/dh); %Indices y, vertical ny2=round(capa2/dh); ny3=round(capa3/dh); ny=ny1+ny2+ny3; %Total de indices y rho=1200;c=1500;k=(c^2)*rho; %Densidad, velocidad y k = bulk modulus dt=dh/c/4; %Diferencia temporal entre pasos %Crea una matriz llena de zeros -> dp las rellana con las densidades rho=zeros(nx,ny); rho(:,ny2+ny3+1:ny)=1100; %Blanda rho(:,ny3+1:ny2+ny3)=1160; %Dura rho(:,1:ny3)=1100; %Blanda rhox=(rho(1:end-1,:)+rho(2:end,:))/2; %Rellenamos bulk x rhoy=(rho(:,1:end-1)+rho(:,2:end))/2; %Rellenamos bulk y imagesc(transpose(rho)); %Representamos rho de la gelatina title([name ': Impedancia acústica']) %Ponemos título pause(4); close; %Matriz presión y velocidades p=zeros(nx,ny); %Presión ux=zeros(nx+1,ny); %Velocidad uy=zeros(nx,ny+1); %Velocidad %Excitación centralfrequency=.5e6; %Frecuencia resonancia transductor a=centralfrequency/(sqrt(pi)/2)*4; %Frecuencia en radianes t=((1:maxtt*maxttt)/(1/dt)-4/a); w=-(exp(-a^2*(t.^2)/2).*(a^2*(t.^2)-1)); h=[]; dh=dh*(100-niveldecompresion)/100; %Comprime %Dibujar por primera vez pcolor((1:nx)*dh,(1:ny)*dh,10*log10(min(max(p'.^2,1d-6),1)) ); shading flat; axis equal; colorbar; title([name ': Time=0 micros']); drawnow %Captura la primera imágen si captureStatus == 1 if captureMovie == 1 120 pause(1); framesCollector(j)= getframe(gcf); %Captura frames close; j= j + 1; else %Se queja por si se quiere ver la pinícula cuando no se ha pedido framesCollector = 'Has dicho que no querías movie!'; end %Se queja por si se quiere ver la pinícula cuando no se ha pedido %Varios util1=round(nx/2-.025/dh):round(nx/2+.025/dh); util2=((1-cos(2*pi*(0:length(util1)-1)./(length(util1)-1)))/2); %Empieza a calcular for ttt=1:maxttt; for tt=1:maxtt; %Cálculo presión p=p-(c.^2)*dt/dh*(diff(ux)+diff(uy')').*rho; %Cálculo velocidad ux(2:nx,:)=ux(2:nx,:)-dt/dh*diff(p)./rhox; uy(:,2: ny)=uy(:,2:ny)-dt/dh*diff(p')'./rhoy; %Condiciones de contorno ux(1,:)=-0*p(1,:)/5000000000; ux(end,:)=0*p(end,:)/5000000000; uy(:,1)=0; uy(:,end)=p(:,end)/5000000000; if nx>length(util1) uy(util1,ny)=uy(util1,ny)+util2.'*w(tt+maxtt*(ttt-1))*dt*10; else uy(:,ny)=uy(:,ny)+w(tt+maxtt*(ttt-1))*dt*10; end %End del for interno %Recording h=[h p(round(nx/2),ny)]; end %End del for externo %Representamos pcolor((1:nx)*dh,(1:ny)*dh,10*log10(min(max(p'.^2,1d-6),1)) ); shading flat; axis equal; colorbar; set(gca,'CLim',[-60 0]); title([name ': Time = ' num2str(round((tt+maxtt*(ttt-1))*1000000*dt)) ' micros']); drawnow %Representa, colorbar, fija ejes, título... %Guarda el frame actual y avanza el puntero una posición if captureMovie == 1 %Solo captura si se pide que se capture captureMovie pause(1); %Esperamos a que se harya refrescado y hayan desaparecido %los efectos especiales de Win 7 framesCollector(j)=getframe(gcf); %Guardamos frame close; %Cierra, molesto pero problemas con los frames j= j + 1; end end return