Este artículo fue publicado como parte del Blogatón de ciencia de datos
La extracción del cráneo es uno de los pasos preliminares en el camino para detectar anomalías en el cerebro. Es el proceso de aislar tejido cerebral de tejido no cerebral a partir de una imagen de resonancia magnética de un cerebro. Esta segmentación del cerebro del cráneo es una tarea tediosa incluso para radiólogos expertos y la precisión de los resultados varía mucho de una persona a otra. Aquí estamos tratando de automatizar el proceso creando una canalización de extremo a extremo donde solo necesitamos ingresar la imagen de resonancia magnética sin procesar y la canalización debe generar una imagen segmentada del cerebro después de realizar el preprocesamiento necesario.
Entonces, ¿qué es una imagen de resonancia magnética?
Para obtener una imagen de resonancia magnética de un paciente, se insertan en un túnel con un campo magnético en su interior. Esto hace que todos los protones del cuerpo se «alineen», por lo que su giro cuántico es el mismo. A continuación, se utiliza un pulso del campo magnético oscilante para interrumpir esta alineación. Cuando los protones vuelven al equilibrio, envían una onda electromagnética. En función del contenido de grasa, la composición química y, lo que es más importante, el tipo de estimulación (es decir, las secuencias) utilizadas para interrumpir los protones, se obtendrán diferentes imágenes. Cuatro secuencias comunes que se obtienen son T1, T1 con contraste (T1C), T2, y INSTINTO.
Desafíos comunes al trabajar con imágenes cerebrales
-
Desafíos de los datos del mundo real
Construir un modelo y lograr una buena precisión en una computadora portátil jupyter es bueno. Pero la mayoría de las veces, un modelo de muy buen rendimiento se comporta muy mal en datos del mundo real. Esto sucede debido a la deriva de datos cuando el modelo ve datos completamente diferentes a los que está entrenado. En nuestro caso, puede suceder debido a diferencias en algunos parámetros o métodos de generación de imágenes de resonancia magnética. Aquí hay un Blog que describe algunos fallos de la IA en el mundo real.
Formulación del problema
La tarea que tenemos aquí es dar una imagen de resonancia magnética en 3D que tenemos para identificar el cerebro y segmentar el tejido cerebral de la imagen completa de un cráneo. Para esta tarea, tendremos una etiqueta de verdad básica y, por lo tanto, será una tarea de segmentación de imágenes supervisada. Nosotros usaremos pérdida de dados como nuestra función de pérdida.
Datos
Echemos un vistazo al conjunto de datos que usaremos para esta tarea. El conjunto de datos se puede descargar desde aquí.
El repositorio contiene datos de 125 participantes, de 21 a 45 años, con una variedad de síntomas psiquiátricos clínicos y subclínicos. Para cada participante, el repositorio contiene:
- Imagen de resonancia magnética ponderada en T1 estructural (sin rostro): esta es la imagen de resonancia magnética ponderada en T1 sin procesar con un solo canal.
- Máscara cerebral: es la máscara de imagen del cerebro o se puede llamar la verdad fundamental. Se obtiene utilizando el método Beast (extracción de cerebro basada en segmentación no local) y aplicando ediciones manuales por expertos en el dominio para eliminar tejido no cerebral.
- Imagen de cráneo despojado: esto se puede considerar como parte del cerebro despojado de la imagen ponderada en T1 anterior. Esto es similar a superponer máscaras a imágenes reales.
La resolución de las imágenes es de 1 mm3 y cada archivo está en formato NiFTI (.nii.gz). Un único punto de datos se parece a esto …
Preprocesando nuestras imágenes Raw
img=nib.load('/content/NFBS_Dataset/A00028185/sub-A00028185_ses-NFB3_T1w.nii.gz')
print('Shape of image=",img.shape)
Imagine imágenes de arriba en 3-D como si tuviéramos 192 imágenes en 2-D de tamaño 256 * 256 apiladas una encima de la otra.
Creemos un marco de datos que contenga la ubicación de las imágenes y las máscaras correspondientes y las imágenes sin calaveras.
#storing the address of 3 types of files
import os
brain_mask=[]
brain=[]
raw=[]
for subdir, dirs, files in os.walk("/content/NFBS_Dataset'):
for file in files:
#print os.path.join(subdir, file)y
filepath = subdir + os.sep + file
if filepath.endswith(".gz"):
if '_brainmask.' in filepath:
brain_mask.append(filepath)
elif '_brain.' in filepath:
brain.append(filepath)
else:
raw.append(filepath)
La señal de campo de polarización es una señal de baja frecuencia y muy suave que corrompe las imágenes de resonancia magnética, especialmente las producidas por resonancias magnéticas antiguas (Imagen de resonancia magnética) máquinas. Los algoritmos de procesamiento de imágenes como la segmentación, el análisis de texturas o la clasificación que utilizan los valores de nivel de gris de los píxeles de la imagen no producirán resultados satisfactorios. Se necesita un paso de preprocesamiento para corregir la señal del campo de polarización antes de enviar imágenes de resonancia magnética corruptas a dichos algoritmos o los algoritmos deben modificarse.
-
Recortar y cambiar el tamaño
Debido a las limitaciones computacionales de ajustar la imagen completa al modelo aquí, decidimos reducir el tamaño de la imagen de resonancia magnética de (256 * 256 * 192) a (96 * 128 * 160). El tamaño del objetivo se elige de tal manera que se captura la mayor parte del cráneo y, después de recortarlo y redimensionarlo, tiene un efecto de centrado en las imágenes.
-
Normalización de la intensidad
La normalización cambia y escala una imagen para que los píxeles de la imagen tengan una media cero y una varianza unitaria. Esto ayuda a que el modelo converja más rápido al eliminar la variación de escala. A continuación se muestra el código para ello.
class preprocessing(): def __init__(self,df): self.data=df self.raw_index=[] self.mask_index=[] def bias_correction(self): !mkdir bias_correction n4 = N4BiasFieldCorrection() n4.inputs.dimension = 3 n4.inputs.shrink_factor = 3 n4.inputs.n_iterations = [20, 10, 10, 5] index_corr=[] for i in tqdm(range(len(self.data))): n4.inputs.input_image = self.data.raw.iloc[i] n4.inputs.output_image="bias_correction/"+str(i)+'.nii.gz' index_corr.append('bias_correction/'+str(i)+'.nii.gz') res = n4.run() index_corr=['bias_correction/'+str(i)+'.nii.gz' for i in range(125)] data['bias_corr']=index_corr print('Bias corrected images stored at : bias_correction/') def resize_crop(self): #Reducing the size of image due to memory constraints !mkdir resized target_shape = np.array((96,128,160)) #reducing size of image from 256*256*192 to 96*128*160 new_resolution = [2,]*3 new_affine = np.zeros((4,4)) new_affine[:3,:3] = np.diag(new_resolution) # putting point 0,0,0 in the middle of the new volume - this could be refined in the future new_affine[:3,3] = target_shape*new_resolution/2.*-1 new_affine[3,3] = 1. raw_index=[] mask_index=[] #resizing both image and mask and storing in folder for i in range(len(data)): downsampled_and_cropped_nii = resample_img(self.data.bias_corr.iloc[i], target_affine=new_affine, target_shape=target_shape, interpolation='nearest') downsampled_and_cropped_nii.to_filename('resized/raw'+str(i)+'.nii.gz') self.raw_index.append('resized/raw'+str(i)+'.nii.gz') downsampled_and_cropped_nii = resample_img(self.data.brain_mask.iloc[i], target_affine=new_affine, target_shape=target_shape, interpolation='nearest') downsampled_and_cropped_nii.to_filename('resized/mask'+str(i)+'.nii.gz') self.mask_index.append('resized/mask'+str(i)+'.nii.gz') return self.raw_index,self.mask_index def intensity_normalization(self): for i in self.raw_index: image = sitk.ReadImage(i) resacleFilter = sitk.RescaleIntensityImageFilter() resacleFilter.SetOutputMaximum(255) resacleFilter.SetOutputMinimum(0) image = resacleFilter.Execute(image) sitk.WriteImage(image,i) print('Normalization done. Images stored at: resized/')
Echemos un vistazo a la arquitectura del modelo.
def data_gen(self,img_list, mask_list, batch_size):
'''Custom data generator to feed image to model'''
c = 0
n = [i for i in range(len(img_list))] #List of training images
random.shuffle(n)
while (True):
img = np.zeros((batch_size, 96, 128, 160,1)).astype('float') #adding extra dimensions as conv3d takes file of size 5
mask = np.zeros((batch_size, 96, 128, 160,1)).astype('float')
for i in range(c, c+batch_size):
train_img = nib.load(img_list[n[i]]).get_data()
train_img=np.expand_dims(train_img,-1)
train_mask = nib.load(mask_list[n[i]]).get_data()
train_mask=np.expand_dims(train_mask,-1)
img[i-c]=train_img
mask[i-c] = train_mask
c+=batch_size
if(c+batch_size>=len(img_list)):
c=0
random.shuffle(n)
yield img,mask
Estamos utilizando una U-Net 3D como nuestra arquitectura. Si ya está familiarizado con la U-Net 2D, esto será muy simple. Primero, tenemos una ruta de contracción a través de un codificador que reduce gradualmente el tamaño de la imagen y aumenta el número de filtros para generar características de cuello de botella. Esto luego se alimenta a un bloque de decodificador que expande gradualmente el tamaño para que finalmente pueda generar una máscara como la salida prevista.
def convolutional_block(input, filters=3, kernel_size=3, batchnorm = True):
'''conv layer followed by batchnormalization'''
x = Conv3D(filters = filters, kernel_size = (kernel_size, kernel_size,kernel_size),
kernel_initializer="he_normal", padding = 'same')(input)
if batchnorm:
x = BatchNormalization()(x)
x = Activation('relu')(x)
x = Conv3D(filters = filters, kernel_size = (kernel_size, kernel_size,kernel_size),
kernel_initializer="he_normal", padding = 'same')(input)
if batchnorm:
x = BatchNormalization()(x)
x = Activation('relu')(x)
return x
def resunet_opt(input_img, filters = 64, dropout = 0.2, batchnorm = True):
"""Residual 3D Unet"""
conv1 = convolutional_block(input_img, filters * 1, kernel_size = 3, batchnorm = batchnorm)
pool1 = MaxPooling3D((2, 2, 2))(conv1)
drop1 = Dropout(dropout)(pool1)
conv2 = convolutional_block(drop1, filters * 2, kernel_size = 3, batchnorm = batchnorm)
pool2 = MaxPooling3D((2, 2, 2))(conv2)
drop2 = Dropout(dropout)(pool2)
conv3 = convolutional_block(drop2, filters * 4, kernel_size = 3, batchnorm = batchnorm)
pool3 = MaxPooling3D((2, 2, 2))(conv3)
drop3 = Dropout(dropout)(pool3)
conv4 = convolutional_block(drop3, filters * 8, kernel_size = 3, batchnorm = batchnorm)
pool4 = MaxPooling3D((2, 2, 2))(conv4)
drop4 = Dropout(dropout)(pool4)
conv5 = convolutional_block(drop4, filters = filters * 16, kernel_size = 3, batchnorm = batchnorm)
conv5 = convolutional_block(conv5, filters = filters * 16, kernel_size = 3, batchnorm = batchnorm)
ups6 = Conv3DTranspose(filters * 8, (3, 3, 3), strides = (2, 2, 2), padding = 'same',activation='relu',kernel_initializer="he_normal")(conv5)
ups6 = concatenate([ups6, conv4])
ups6 = Dropout(dropout)(ups6)
conv6 = convolutional_block(ups6, filters * 8, kernel_size = 3, batchnorm = batchnorm)
ups7 = Conv3DTranspose(filters * 4, (3, 3, 3), strides = (2, 2, 2), padding = 'same',activation='relu',kernel_initializer="he_normal")(conv6)
ups7 = concatenate([ups7, conv3])
ups7 = Dropout(dropout)(ups7)
conv7 = convolutional_block(ups7, filters * 4, kernel_size = 3, batchnorm = batchnorm)
ups8 = Conv3DTranspose(filters * 2, (3, 3, 3), strides = (2, 2, 2), padding = 'same',activation='relu',kernel_initializer="he_normal")(conv7)
ups8 = concatenate([ups8, conv2])
ups8 = Dropout(dropout)(ups8)
conv8 = convolutional_block(ups8, filters * 2, kernel_size = 3, batchnorm = batchnorm)
ups9 = Conv3DTranspose(filters * 1, (3, 3, 3), strides = (2, 2, 2), padding = 'same',activation='relu',kernel_initializer="he_normal")(conv8)
ups9 = concatenate([ups9, conv1])
ups9 = Dropout(dropout)(ups9)
conv9 = convolutional_block(ups9, filters * 1, kernel_size = 3, batchnorm = batchnorm)
outputs = Conv3D(1, (1, 1, 2), activation='sigmoid',padding='same')(conv9)
model = Model(inputs=[input_img], outputs=[outputs])
return model
Luego entrenamos el modelo usando el optimizador de Adam y la pérdida de dados como nuestra función de pérdida …
def training(self,epochs):
im_height=96
im_width=128
img_depth=160
epochs=60
train_gen = data_gen(self.X_train,self.y_train, batch_size = 4)
val_gen = data_gen(self.X_test,self.y_test, batch_size = 4)
channels=1
input_img = Input((im_height, im_width,img_depth,channels), name="img")
self.model = resunet_opt(input_img, filters=16, dropout=0.05, batchnorm=True)
self.model.summary()
self.model.compile(optimizer=Adam(lr=1e-1),loss=focal_loss,metrics=[iou_score,'accuracy'])
#fitting the model
callbacks=callbacks = [
ModelCheckpoint('best_model.h5', verbose=1, save_best_only=True, save_weights_only=False)]
result=self.model.fit(train_gen,steps_per_epoch=16,epochs=epochs,validation_data=val_gen,validation_steps=16,initial_epoch=0,callbacks=callbacks)
Después de entrenar durante 60 épocas, obtuvimos un iou_score de validación de 0.86.
Echemos un vistazo a cómo funcionó nuestro modelo. Nuestro modelo predecirá simplemente la máscara. Para obtener la imagen sin calavera, necesitamos superponerla en la imagen Raw para obtener la imagen sin calavera …
Mirando las predicciones podemos decir que aunque es capaz de identificar el cerebro y segmentarlo, no se acerca a la perfección. En este punto, podemos sentarnos con un experto en el dominio para identificar qué pasos adicionales de preprocesamiento se pueden realizar para mejorar la precisión. Pero en cuanto a este post, lo concluiré aquí. Por favor, siga enlace1 y / o enlace2 si quieres saber mas …
Conclusión:
Me alegro de que hayas llegado al final. Espero que esto te ayude a comenzar con la segmentación de imágenes en imágenes 3D. Puede encontrar el enlace de Google Colab que contiene el código. aquí. No dude en agregar sugerencias o consultas en la sección de comentarios. ¡Que tenga un lindo día!
Los medios que se muestran en este artículo sobre la eliminación del cráneo no son propiedad de DataPeaker y se utilizan a discreción del autor.
