1/4

I. Introducción

Las imágenes de ultrasonido se utilizan como la principal forma de imagenología en diversas especialidades médicas, en su mayoría para diagnosticar padecimientos en tejidos blandos como riñones o hígado. También es usado para dar seguimiento al desarrollo fetal durante el embarazo. Recientemente, se han utilizado para visualizar el tejido intraarticular, el cartílago, los ligamentos; y para diagnosticar la fractura de huesos que, en algunos casos puede sustituir a los tradicionales rayos X. [1] Otra alternativa viable, es la navegación y guía transoperatoria en cirugías ortopédicas asistidas por computadora (CAOS por sus siglas en inglés), donde las imágenes de fluoroscopía o de tomografía computarizada son las más usadas, esto debido a los beneficios que el ultrasonido brinda como: la portabilidad, el bajo costo, así como la posibilidad de adquirir imágenes libres de radiaciones.

En muchos procedimientos quirúrgicos, no basta con hacer una visualización de las imágenes de ultrasonido, es necesario contar con modelos anatómicos virtuales, para realizar una segmentación de los órganos de interés. En particular para este estudio serán los huesos, específicamente la meseta tibial.

La segmentación del hueso, y en general la segmentación de estructuras a partir de imágenes de ultrasonido, es una tarea difícil dadas las características físicas de este, principalmente por la baja relación entre señal y ruido. En particular, las estructuras hiperecoicas que poseen los huesos (i. e. con alta densidad, que produce ecos ultrasónicos de alta amplitud), producen una gran sombra acústica, reverberaciones, y una respuesta distinta dependiendo del ángulo de incidencia de las ondas acústicas.

En este artículo reportamos la segmentación de la superficie ósea de la meseta tibial. Se implementó el aprendizaje de máquina basado en un clasificador supervisado de Bayes, con el objetivo de demostrar la viabilidad del uso de esta metodología en la segmentación de volúmenes de ultrasonido, los cuales se obtuvieron mediante la técnica conocida como “ultrasonido 3D a manos libres” descrita en Image tracking and volume reconstruction of medical ultrasound. [2] El objetivo final es poder realizar el registro transoperatorio en una cirugía ortopédica asistida por computadora, con el modelo virtual y el plan preoperatorio del paciente, como se reportó en Improved surface-based registration of CT and intraoperative 3D ultrasound of bones. [3]

2/5

II. Detección de la superficie ósea en ultrasonido 3D

Las estructuras óseas aparecen en las imágenes de ultrasonido como regiones con bordes brillantes seguidas de una región obscura, resultado de la sombra acústica generada por el hueso como se observa en la figura 1. La alta intensidad del reflejo acústico provocado por el hueso se observa en las imágenes como una línea de espesor entre 2 y 4mm. Característica que debe ser tomada en cuenta para que las dimensiones de la superficie extraída correspondan con el tamaño original del hueso del paciente. Jain et al., [4] demuestran que la cortical del hueso se encuentra en el máximo de intensidad del gradiente generado por la respuesta del hueso al ultrasonido. El resultado activa un detector que realza la línea media para que tenga una representación más precisa de la cortical ósea en imágenes de ultrasonido que los métodos basados en detección de bordes.

Nuestro método de segmentación toma en cuenta alteraciones como el ruido speckle y las reverberaciones, para minimizarlos se desarrollaron dos procesos: el realce de estructuras óseas y la posterior segmentación sobre la imagen realzada, como se describe a continuación.

II.1. Realce de superficies en ultrasonido 3D

La ubicación de una cresta local en una función se da cuando la segunda derivada es igual a cero. Sea I un volumen de ultrasonido, este se puede expresar como una función en tres dimensiones multivariable, por lo que las segundas derivadas direccionales en x, y y z del volumen describen la variación del gradiente de intensidad de los voxeles localmente. La segunda derivada direccional en un punto x̂ =(x,y,z) del volumen se representa como matriz hessiana (Ec. 1, apéndice 1), donde cada elemento representa la derivada parcial de segundo orden de la imagen en cada dirección (x,y,z). Una propiedad de esta matriz es que es simétrica con eigenvalores reales y eigenvectores ortogonales, por lo tanto, invariante a rotaciones.

La matriz hessiana se construye a partir del cálculo de la segunda derivada parcial direccional, que puede ser aproximada mediante la convolución del volumen I con la segunda derivada de la función Gaussiana normalizada G (Ec. 2 y 3, apéndice 1), donde la desviación estándar σ de G es usado como un factor de escala, el cual se asignará dependiendo del tamaño en la imagen de las estructuras que se quieran realzar.

La segunda derivada de la imagen aplicada sobre un voxel, proporciona información del tipo de estructura a la que pertenece, pueden ser formas tubulares, esféricas o de superficie. Esta información se obtiene del eigenanálisis de la matriz hessiana generada para cada voxel en el volumen (Tabla 1, apéndice 1). Se eligen solo los voxeles que cumplan con la condición de pertenecer a una superficie, un valor de “superficialidad” les es asignado (Ec. 4, apéndice 1) y como resultado se obtiene una imagen 3D donde cada voxel representa una medida de “superficialidad” para una escala σ dada, la cual llamaremos imagen realzada.

II.2. Clasificación de las superficies de los huesos

Posteriormente a partir de la imagen realzada, la segmentación de la cortical del hueso se hace usando un clasificador supervisado de Bayes, el cual asigna cada voxel en una de las posibles clases del conjunto C=c₁,c₂,···c_m definidas previamante. La asignación de un voxel a una clase se da cuando la probabilidad de pertenecer a la clase C_k dado un vector que caracteriza a este, (X̂) es máxima sobre las demás clases. Esta probabilidad se obtiene a partir del teorema de Bayes definido en la ecuación 6 en el apéndice 1.

Si asumimos que la distribución N-dimensional para cada una de las clases es Gaussiana, se puede obtener el discriminante de Bayes (Ec. 7, apéndice 1) para cada una de las clases k, el cual representará la probabilidad del vector X̂ de pertenencia a la clase k.

Para las imágenes de ultrasonido de hueso se definieron tres clases, que corresponden a diferentes regiones y claramente distintas entre sí: el tejido, que es todo lo visible por encima del hueso; la cortical ósea que se observa como la línea más brillante en la imagen; y la sombra acústica, que es la sección obscura por debajo del hueso. El vector de características para cada voxel se definió como: la intensidad del voxel en el volumen realzado (E(x,y,z)), y los tres primeros momentos estadísticos del volumen realzado (M_e1 (x,y,z), M_e2 (x,y,z), M_e3 (x,y,z)) definidos en las ecuaciones 8, 9 y 10 respectivamente, incluidas en el apéndice 1.

Una vez construido el vector de características para cada voxel en el volumen se procede a entrenar el clasificador supervisado, para esto se necesita un volumen de entrenamiento que tenga tres conjuntos de voxeles etiquetados por un radiólogo experto con la clase a la que pertenecen. Con este volumen se encuentran la media (μ_k) y la matriz de covarianza (Σ_k) para cada clase k.

El entrenamiento del sistema solo se realiza una vez para imágenes pertenecientes a un mismo protocolo de adquisición y se puede realizar fuera de línea. En la etapa de clasificación, se construye una imagen de probabilidad I_p (x,y,z), donde el valor de cada voxel corresponde a la probabilidad de pertenecer a la clase hueso. Si el valor de probabilidad es mayor para otra clase se le asignará a dicho voxel un valor de cero. Finalmente, a partir de la imagen de probabilidad se selecciona un porcentaje mínimo de voxeles con el valor de probabilidad más alto de pertenecer a la cortical ósea, a este conjunto se designará como el conjunto de semillas, las cuales son usadas para inicializar un algoritmo de crecimiento de regiones, [5] a lo largo de todos los voxeles vecinos cuyo valor de probabilidad es I_p>∈. Donde ∈ es un valor que se asigna de forma experimental con el cual podemos modular que tan “grueso” se quiere la sección de la cortical del hueso. En esta parte, el algoritmo de crecimiento de regiones se usa para evitar que artefactos como las reverberaciones sean clasificadas como hueso, ya que estas presentan una característica similar.

II. 3. Adquisición de imágenes de ultrasonido 3D de la tibia

Para validar el método de segmentación se utilizó un maniquí realista, que corresponde a la tibia como se muestra en la figura 6. Se adquirieron 12 volúmenes de ultrasonido de distintas partes que involucran a la meseta tibial. Las imágenes de ultrasonido del maniquí se adquirieron con la técnica de ultrasonido navegado 3D de mano libre descrita en Torres et al., [2] En la figura 2 se observan tres planos ortogonales de ultrasonido del maniquí de la tibia.

Se usó un modelo virtual de imágenes de tomografía computarizada del maniquí, para evaluar los resultados de la segmentación de la superficie de la tibia en el volumen de ultrasonido, como se describe a continuación.

II. 4. Modelo virtual de tomografía computarizada (TC) de la tibia

Uno de los métodos de tomografía más utilizados son los métodos de tomografía computarizada por rayos X (TC) como lo señalan Fishman & Jefrey. [6] La TC funciona obteniendo imágenes de secciones del objeto en forma secuencial. Existen diversas configuraciones para la adquisición de estas secciones, por un lado, se cuenta con la configuración helicoidal utilizada comúnmente en la medicina, como se muestra en la figura 3, y la utilizada en la industria, la cual utiliza una mesa giratoria como el sistema de adquisición en los sistemas de microtomografía, figura 4. [7]

Los TC proveen un grupo de imágenes I, donde cada una representa una sección 2D del volumen, el grupo de todas define el volumen digitalizado V ∈ Z³. Al igual que con las imágenes donde el concepto de píxel[8] se utiliza, para definir la menor unidad homogénea que contiene información de la imagen. Los volúmenes V estan formados por voxeles ϒ ∈ V y son los elementos procesables más pequeños. La resolución espacial de nuestro volumen depende de la resolución del sensor y la distancia recorrida entre cada toma, es común tener diferentes resoluciones en cada dirección espacial.

Los valores guardados en cada voxel y el proceso en el que la imagen tomográfica se forma pueden entenderse como un símil a la transparencia y translucencia de objetos con respecto a la luz visible, donde los materiales presentan diferentes valores de radiodensidad que indican que tan opaco es el material a los rayos X como lo señalan Novelline & Squiere. [9]

La profundidad de bits [10] del sensor y los materiales del objeto a adquirir definen el rango dinámico del volumen, cada sensor tiene distintas capacidades para representar varias radiodensidades de forma discreta. Por ejemplo, un sensor de 8 Bits puede guardar 256 diferentes radiodensidades como valores de gris. Y la complejidad de los materiales a adquirir en una toma define que porciones del rango dinámico son usadas para representar cada material. En la figura 5 se muestra una imagen obtenida por técnicas de CT con su rango dinámico. [10]

Para construir el modelo virtual de la tibia se usó un modelo de CT basado en un maniquí de tibia. Este fue construido con hueso tibial sintético (ERP #1117-42, Sawbones Inc., Vashon, WA, lE.E.U.U.) sumergido en un hidrogel hecho del alcohol de polivinilo (PVA), diluido en el 95% de agua. En la figura 6 se muestra una imagen del hueso sintético y del maniquí final cubierto de hidrogel.

El maniquí fue digitalizado con un MicroCT, modelo Nikon Metrology XTH 225, con resolución de imagen de 2048 x 2048. Se configuró el sistema, para obtener imágenes utilizando 220 kV y 61 µA, con un numero de proyecciones de 3142. Esta configuración permitió obtener un volumen final compuesto por 1042 × 1250 × 3201 voxeles donde cada voxel tiene una resolución de 0.115 mm³.

La superficie ósea se obtuvo etiquetando el voxel que contiene información del hueso sintético con el valor de 1 y el resto de los voxeles se les asignó un valor de 0. El proceso de etiquetado se realizó utilizando un algoritmo de agrupamiento k-means. [11] El volumen del hueso obtenido del CT se muestra en la figura 7.

3/5

III. Resultados

Para cada volumen segmentado a partir de las imágenes de ultrasonido se montó el modelo virtual de CT y se obtuvo la distancia mínima al volumen de referencia para cada uno de los voxeles segmentados. En promedio, se registró una distancia máxima de 2.031 mm, sin embargo, el promedio de la distancia media es 0.21 mm con una desviación estándar de 0.17 mm. Esto muestra que muy pocos de los voxeles segmentados están fuera de la cortical del hueso. En la figura 8 (sección superior) se muestra la imagen original de US, la segmentación de referencia en amarillo y la segmentación del volumen de US sobrepuesta en azul.

Una vez segmentado el volumen, y a partir de este, se construye una malla la cual representa la anatomía del paciente en la mesa de operaciones. El modelo preoperatorio se pone en el mismo sistema de referencias que la segmentación para poder inicializar el proceso de registro, con el algoritmo Iterative Closest Points. [12] La figura 8 (sección inferior) muestra el resultado del registro de las dos mallas, la segmentada en azul y la del plan preoperatorio en amarillo.

4/5

IV. Conclusiones

En este trabajo se presenta un método, para la detección automática de la superficie de los huesos en imágenes de ultrasonido 3D, basado en la detección de estructuras 3D: superficies, líneas, esferas y en la clasificación automática de la superficie de los huesos utilizando un clasificador de Bayes. La validación del método de segmentación en el maniquí de tibia produjo resultados con errores de exactitud suficientemente pequeños como para implementar este método en la detección de huesos con fines quirúrgicos, aunque una validación extendida a otros huesos será necesaria.

Anexo 1

Configuración de la matriz hessiana:

cada elemento I_ij representa la derivada parcial de segundo orden de la imagen en la dirección ij en la posición x̂, esto es, I_xx=∂² I / ∂ x², I_xy=∂² I / ∂ x ∂y,⋅⋅⋅ y así sucesivamente. Cuando i=j entonces I_ij=I_ji, esta propiedad hace que la matriz hessiana sea simétrica con eigenvalores reales y eigenvectores ortogonales, por lo tanto, invariante a rotaciones. Las segundas derivadas se aproximan con la convolución del volumen I con la segunda derivada de la función Gaussiana normalizada G, de esta manera para la dirección ij, el elemento I_ij de la matriz hessiana se define como:

donde:

con σ el valor de la desviación estándar. El valor de σ es también utilizado como un factor de escala, el cual se asignará dependiendo del tamaño en la imagen de las estructuras que se quieran realzar.

El eigenanálisis para cada voxel en el volumen resulta en la obtención de tres eigenvalores λ₁, λ₂ y λ₃. Sin pérdida de generalidad, sean λ₁ ≥ λ₂ ≥ λ₃, se puede encontrar un conjunto de condiciones discriminantes de pertenencia a un tipo de estructura, la cual puede ser tubular, esférica o de superficie. [6] Estas condiciones se muestran a continuación en la tabla 1:

La medida de superficialidad para cada voxel se puede asignar a partir de la siguiente ecuación:

donde:

los parámetros α y γ fueron escogidos experimentalmente, obteniendo α = 0.5 y γ = 1.

Con la imagen de probabilidades y basados en el teorema de Bayes se puede obtener la probabilidad de pertenecer a una clase k dado un vector de características:

donde:

P(C_k) es la probabilidad a priori de ocurrencia para la clase k.

P(X̂│C_k) es la probabilidad condicional de X̂ dado C_k.

X̂ el vector û = u₁, u₂,⋅⋅⋅u_n que caracteriza un voxel del volumen.

Si asumimos que la distribución N-dimensional para cada una de las clases es Gaussiana, el discriminante de Bayes para cada una de las clases k se define como:

donde:

A es la dimensión del vector X̂

μ_k es un vector N-dimensional que representa la media de cada una de las características del conjunto de entrenamiento.

Σ_k es una matriz de N × N la matriz de covarianza para C_k.

|Σ_k| el determinante de la matriz de covarianza.

Los tres primeros momentos estadísticos se definen como:

El cálculo de los momentos estadísticos se realizó dentro de una ventana de tamaño fijo escogido como tres veces el tamaño ocupado por la respuesta del hueso en el ultrasonido, esto con el fin de poder obtener la relación entre los voxeles pertenecientes a la cortical del hueso y los voxeles vecinos a esta.

5/5

Bibliografía

[1] A. L. Waterbrook, S. Adhikari, U. Stolz, et al., “The accuracy of point-of-care ultrasound to diagnose long bone fractures in the ED.” en The American journal of emergency medicine, vol. 31, no. 9, pp. 1352-1356. 2013.

[2] F. Torres, Z. Fanti, E. Lira, et al., “Image tracking and volume reconstruction of medical ultrasound,” en Revista Mexicana de Ingenieria Biomedica, vol. 33, no. 2, pp. 101-115, 2012.

[3] Z. Fanti, E. Torres, A. Hazan-Lasri, et al,. “Improved surface-based registration of CT and intraoperative 3D ultrasound of bones,” en Journal of healthcare engineering, 2018.

[4] A. K. Jain, & R. H. Taylor, “Understanding bone responses in B-mode ultrasound images and automatic bone surface extraction using a Bayesian probabilistic framework,” in Medical imaging 2004: ultrasonic imaging and signal processing, vol. 5373, pp. 131-142. International Society for Optics and Photonics. 2004, April.

[5] R., Adams, & L. Bischof, “Seeded region growing,” en IEEE Transactions on pattern analysis and machine intelligence, vol.16, no. 6, pp. 641-647, 1994.

[6] E. K. Fishman, & R. B. Jeffrey, Spiral CT: Principles, Techniques, and Clinical Applications. Raven Press, 1995.

[7] E. L. Ritman, “Current Status of Developments and Applications of Micro-CT,” en Annual Review of Biomedical Engineering, vol. 13, no. 1, pp.531–552, 2011.

[8] J. D. Foley, & A. Van Dam, Fundamentals of interactive computer graphics. Reading, Mass, Addison-Wesley Pub. Co, 1983.

[9] R. A. Novelline, & L. F. Squire, Squire’s Fundamentals of Radiology. La Editorial, UPR, 2004.

[10] “An In-Depth Look at Bit Depth," Blog Post Teledyne Lumenera. (n.d.). [En línea]. Disponible en https://www.lumenera.com/blog/bit-depth [Consultado en august 9, 2021].

[11] G. A. F. Seber, Multivariate Observations. John Wiley & Sons, 2004.

[12] P. J. Besl, & N. D. McKay, “Method for registration of 3-D shapes. In Sensor fusion IV: control paradigms and data structures,” en International Society for Optics and Photonics. vol. 1611, pp. 586-606, 1992, April.