DETECCIÓN DE NEUMONÍA EN RADIOGRAFÍAS: 2 – TRANSFORMACIÓN Y SEGMENTACIÓN

Publicado por @JoseMa_deCuenca el March 8, 2020, 9:19 p.m.
En esta entrada vamos a ver como procesar imágenes de forma efectiva, para concentrar los datos y reducir los recursos necesarios para el aprendizaje de los algoritmos que usaremos posteriormente. La **transformación** de imágenes digitales consiste en su remuestreado, es decir la introducción de cambios en su resolución (número de píxeles manteniendo o no su relación de aspecto) y/o en su número de canales (colores). Existen utilidades de código abierto para el tratamiento masivo de imágenes, como por ejemplo imagemagick (imagemagick.org), que permiten procesar de forma efectiva carpetas enteras, homogeneizando su resolución, unificando su número de canales, adoptando una única relación de aspecto, etc. Pero aunque probablemente aceleren los tiempos de proceso, se ha preferido utilizar las librerías de python para ello, para obtener mayor control de las transformaciones realizadas y también para poder realizar pruebas de rendimiento y optimizar el proceso. La **segmentación** es un proceso de preprocesado de imágenes que permite reducir la cantidad de datos manteniendo aquella información relevante para nosotros. Se basa en la agrupación de los píxeles que tienen un valor de intensidad similar para separar regiones completas de la imagen, ocultando las superficies no deseadas con una máscara. Se aplica al reconocimiento de objetos, su oclusión, detección de límites, clasificación de imágenes o compresión. La segmentación de imágenes se basa en dos propiedades de los valores de intensidad: la discontinuidad y la similitud. La discontinuidad hace una partición de la imagen basada en los cambios en el valor de intensidad nítida, mientras que la propiedad de similitud hace una partición de la imagen en regiones que son similares a los criterios especificados. En este caso se usa para reconocer el contorno de clasificación y suprimir el resto de la imagen. Para ello se utiliza una red neuronal, pero en vez de construirla de cero se recurre a la transferencia de aprendizaje desde una red de arquitectura U-net ya pre entrenada para este uso: el reconocimiento de pulmones en radiografías. Esta red ha sido originalmente entrenada con contornos pulmonares dibujados sobre imágenes reales. Con los pesos utilizados es capaz de reconocer estos órganos a pesar de las diferencias anatómicas de los sujetos, aunque se observa un menor rendimiento en individuos enfermos. No obstante, debemos modificar el código para que ahora oculte toda la imagen menos los pulmones, sin marcarlos. Comenzamos definiendo las funciones a utilizar # SEGMENTADOR PULMONAR CON TRANSFERENCIA DE APRENDIZAJE DE RED NEURONAL PREEXISTENTE EN LA RED: # https://github.com/imlab-uiip/lung-segmentation-2d/blob/master/Demo/demo.py # El código incluido a continuación está inspirado en la demo del repositorio # Importación de librerías y carga de funciones import numpy as np import pandas as pd from keras.models import load_model from keras.preprocessing.image import ImageDataGenerator from skimage import morphology, io, color, exposure, img_as_float, transform from matplotlib import pyplot as plt import os from PIL import Image import time import cv2 from skimage.transform import resize # Funcion para cargar y pre procesar datos. Usa formato im_shape de 256x256, necesario para RN # Si se cargan las imagenes con 256x256 para segmentar def loadDataGeneral(df, path, im_shape): X = [] # matriz para las imagenes for i, item in df.iterrows(): # bucle de carga img = img_as_float(io.imread(item[0])) img = transform.resize(img, im_shape) img = exposure.equalize_hist(img) img = np.expand_dims(img, -1) # expande la forma de la matriz X.append(img) X = np.array(X) # Normalizacion de las imagenes de carga (no de las mascaras) X -= X.mean() X /= X.std() return X # Devuelve array de imagenes a procesar # Funcion para componer imagen coloreando cada pulmon en la imagen con la mascara predicha def mascara(img, mask): rows, cols = img.shape color_mask = np.zeros((rows, cols, 3)) color_mask[mask == 1] = [0, 0, 1] # fondo de color azul img_color = np.dstack((img, img, img)) img_hsv = color.rgb2hsv(img_color) color_mask_hsv = color.rgb2hsv(color_mask) img_hsv[..., 0] = color_mask_hsv[..., 0] img_hsv[..., 1] = color_mask_hsv[..., 1] img_mascara = color.hsv2rgb(img_hsv) return img_mascara # Retorna una imagen con la mascara Las imágenes transformadas y enmascaradas se almacenarán sobre una nueva estructura de carpetas en disco, con su nombre original, replicando la original del dataset utilizado. De esta manera solamente será necesario ejecutar una vez este código, que es bastante pesado, aunque no tanto como para no poder abordarlo de forma doméstica. Tarda unas cuatro horas sobre un i7 de octava generación, y eso a pesar de usar transferencia de aprendizaje. Como las imágenes enmascaradas mantienen su nombre, únicamente será necesario cambiar la ruta de la primera carpeta para utilizar el nuevo dataset segmentado con el resto de cuadernos, reduciendo la carga de programación, además de las veces que debe procesarse. El código trabaja por lotes, uno para cada carpeta de datos del datset original porque se ha comprobado que el procesado resulta más rápido que atacando todo el conjunto a la vez. Esto probablemente es debido a la memoria disponible, y puede que sea necesario usar lotes más pequeños en caso de contar con pocos recursos. Estos lotes se generan a partir del dataframe de características que se generó para la exploración inicial de los datos. Por tanto, tras la definir las funciones a utilizar, en el nuevo cuaderno para datset segmentado se genera el dataframe con el extractor de características que hemos recogido en el primer código, aunque también podría cargarse desde el archivo csv en el que se guardó. A continuación partiremos este dataframe en varios archivos csv, filtrando por conjuntos (entrenamiento, prueba, validación) y por tipo de diagnóstico (normal o pneumonía), de manera que coincidan con las rutas de las carpetas donde se guardan las imágenes del dataset original. También se replica la estructura de carpetas en caso de que no exista, a partir del directorio del usuario, pero colgando de una carpeta chest_xray_seg, para que una vez almacenadas sobre ella, se diferencien de las no segmentadas. # Genero ficheros CSVs separados para carpeta de datos que replicara el segmentador # Divido los conjuntos para procesar por partes ## CONJUNTO DE VALIDACION VAL_NOR_file = (df_file[(df_file['conjunto']=='VAL') & (df_file['diagnostico']=='NORMAL')]['fichero'].iloc[:]) VAL_NOR_file.to_csv ('./data/chest_xray/VAL_NOR_file.csv', index = None, header=True) VAL_PNE_file = (df_file[(df_file['conjunto']=='VAL') & (df_file['diagnostico']=='PNEUMONIA')]['fichero'].iloc[:]) VAL_PNE_file.to_csv ('./data/chest_xray/VAL_PNE_file.csv', index = None, header=True) ## CONJUNTO DE PRUEBA TEST_NOR_file = (df_file[(df_file['conjunto']=='TEST') & (df_file['diagnostico']=='NORMAL')]['fichero'].iloc[:]) TEST_NOR_file.to_csv ('./data/chest_xray/TEST_NOR_file.csv', index = None, header=True) TEST_PNE_file = (df_file[(df_file['conjunto']=='TEST') & (df_file['diagnostico']=='PNEUMONIA')]['fichero'].iloc[:]) TEST_PNE_file.to_csv ('./data/chest_xray/TEST_PNE_file.csv', index = None, header=True) ## CONJUNTO DE ENTRENAMIENTO TRAIN_NOR_file = (df_file[(df_file['conjunto']=='TRAIN') & (df_file['diagnostico']=='NORMAL')]['fichero'].iloc[:]) TRAIN_NOR_file.to_csv ('./data/chest_xray/TRAIN_NOR_file.csv', index = None, header=True) # Añado condicion especial de imagenes en escala de grises porque esta red no permite RGB TRAIN_PNE_file = (df_file[(df_file['conjunto']=='TRAIN') & (df_file['diagnostico']=='PNEUMONIA') & (df_file['color']=='L')]['fichero'].iloc[:]) TRAIN_PNE_file.to_csv ('./data/chest_xray/TRAIN_PNE_file.csv', index = None, header=True) # Creo estructura de carpetas si no existen ya try: os.stat('./data/chest_xray_seg') except: os.mkdir('./data/chest_xray_seg') try: os.stat('./data/chest_xray_seg/test') except: os.mkdir('./data/chest_xray_seg/test') try: os.stat('./data/chest_xray_seg/test/NORMAL') except: os.mkdir('./data/chest_xray_seg/test/NORMAL') try: os.stat('./data/chest_xray_seg/test/PNEUMONIA') except: os.mkdir('./data/chest_xray_seg/test/PNEUMONIA') try: os.stat('./data/chest_xray_seg/train') except: os.mkdir('./data/chest_xray_seg/train') try: os.stat('./data/chest_xray_seg/train/NORMAL') except: os.mkdir('./data/chest_xray_seg/train/NORMAL') try: os.stat('./data/chest_xray_seg/train/PNEUMONIA') except: os.mkdir('./data/chest_xray_seg/train/PNEUMONIA') try: os.stat('./data/chest_xray_seg/val') except: os.mkdir('./data/chest_xray_seg/val') try: os.stat('./data/chest_xray_seg/val/NORMAL') except: os.mkdir('./data/chest_xray_seg/val/NORMAL') try: os.stat('./data/chest_xray_seg/val/PNEUMONIA') except: os.mkdir('./data/chest_xray_seg/val/PNEUMONIA') ### GENERADOR IMAGENES SEGMENTADAS PARA CADA ARCHIVO # Defino ubicacion para el archivo de rutas a procesar, a partir del raiz del usuario csv_path = './data/chest_xray/VAL_NOR_file.csv' path = '.' df = pd.read_csv(csv_path) # dataframe con las rutas del fichero ## Define el formato de carga de las imagenes usando la funcion anterior im_shape = (256, 256) # Respeto el original de la RN de segmentación X = loadDataGeneral(df, path, im_shape) # arrays de numpy cargados a partir de las rutas con el formato 256x256 n_test = X.shape[0] # Toma el numero de datos a tratar inp_shape = X[0].shape # Averigua las dimensiones de las imagenes segun la primera (256,256,1) ### MODELO DE SEGMENTACION n_test = X.shape[0] # Toma el numero de datos a tratar en las carpetas de la ruta # Load model - Transferencia de aprendizaje model_name = './checkpoints/Segmentation_model.hdf5' # modelo entrenado UNet = load_model(model_name) # funcion de keras para cargar modelo # Defino generador para extraer secuencia de pixeles de las imagenes test_gen = ImageDataGenerator(rescale=1.) # Indica reescalado 1. i = 0 # Hace predicciones de mascara a partir del flujo de pixeles alimentado por el generador para cada imagen for xx in test_gen.flow(X, batch_size=1): img = exposure.rescale_intensity(np.squeeze(xx), out_range=(0,1)) # recompone una imagen a partir de sus plixeles, reescalando su intensidad pred = UNet.predict(xx)[..., 0].reshape(inp_shape[:2]) # Usa el modelo cargado para predecir la mascara de la imagen # Crea una mascara binaria a partir de la predicccion pr = pred < 0.5 # Si la prediccion es mayor de 0.5 almacena True # Elimino las regiones pequeñas en las imagenes pr = morphology.remove_small_objects(pr, 0.02 * np.prod(im_shape)) pr = morphology.remove_small_holes(pr, 0.02 * np.prod(im_shape)) # Compongo imagenes enmascaradas if i < n_test: # plt.imshow(pred) # Dibuja la mascara predicha # plt.imshow(mascara(img, pr, pr)) # Dibuja la imagen aplicando la mascara predicha plt.imshow(mascara(img, pr)) # Dibuja la imagen aplicando la mascara predicha # Elimino marcas de ejes plt.xticks([]) plt.yticks([]) # Genero ruta para guardar la nueva imagen ruta = str(df.iloc[i,0]) img_seg = ruta.replace('chest_xray', 'chest_xray_seg') # Guardo la nueva imagen plt.savefig(img_seg) # Graba la figura compuesta de las dos print(img_seg) i += 1 if i == n_test: break # plt.show() # No represento cada una, con la ultima es suficiente Aunque los tiempos de procesado obtenidos con este código son elevados (algo más de cuatro horas), debido al alto número de imágenes, merece la pena porque solo se ejecutará una vez, y se logra una gran reducción del volumen de datos (unas 10 veces, de 1,15 Gb a 111Mb Como se evidenció tras el análisis exploratorio, el conjunto de entrada tiene un sesgo en el número de datos de cada clase que puede trasladarse a los algoritmos desvirtuando sus resultados de rendimiento. Por ello se han realizado pruebas introduciendo un remuestreo simple de los datos de origen, basado en submuestreo, de manera que para el conjunto de entrenamiento se utilice idéntico número de imágenes de una y otra clase. El resampling mantiene los nombres de variables que se usan para alimentar los algoritmos, de manera que el resto de código pueda ser reutilizado. ## RESAMPLING DE LOS DATOS # Reduzco el dataset df_file para eliminar el sesgo en los datos # Homogenizo DIAGNOSTICO para que las clases NORMAL y PNEUMONIA tengan igual numero # Creo subconjuntos de entrenamiento filtrando dataset original por las clases de la caracteristica df_TR_NOR = df_file[(df_file['diagnostico'] == 'NORMAL') & (df_file['conjunto']=='TRAIN')] df_TR_PNE = df_file[(df_file['diagnostico'] == 'PNEUMONIA') & (df_file['conjunto']=='TRAIN')] # Dejo el resto de subconjuntos del dataset en otro dataframe df_RESTO = df_file[df_file['conjunto'] !='TRAIN'] # Cuento datos de cada clases count_NOR = df_TR_NOR.shape[0] count_PNE = df_TR_PNE.shape[0] print('Datos por clase de diagnostico en entrenamiento:') print('count_NOR: ', count_NOR) print('count_PNE: ', count_PNE) # Tomo solo el mismo numero de elementos de la clase mas numerosa que de la clase minoritaria df_TR_PNE_cut = df_TR_PNE.sample(count_NOR) # Recompongo el dataset df_file = pd.concat([df_TR_NOR, df_TR_PNE_cut, df_RESTO], axis=0) # Reviso analisis estadístico con Tabla cruzada de casos print('\nDatos por clase y conjunto:') tab = pd.crosstab(df_file['diagnostico'], df_file['conjunto']) print(tab) Aunque no se ha utilizado, existe una libreria Python (imblearn) específica para el remuestreo equilibrado de un conjunto de datos, compatible con scikit-learn, que permite agrupar los registros de la clase mayoritaria antes de realizar el submuestreo para preservar la información al alejarse de la zona de frontera. También permite introducir ruido o pequeñas variaciones en los duplicados de los datos de la clase minoritaria, para reducir el sobreajuste, aunque no sus costes de proceso. La otra técnica de remuestreo que hubiera sido posible aplicar es un ensemble de diversos conjuntos tomados aleatoriamente, de manera que los modelos entrenados sobre cada uno de ellos proporcione un clasificador no sesgado diferente. Su combinación sería un clasificador general ideal, aplicable al conjunto completo. En el caso de los clasificadores basados en redes neuronales, inestables por naturaleza, las propias diferencias que se producirán en los modelos hacen de este un método apropiado. Las diferentes opciones para realizar los subconjuntos de datos dependen de si se permitirá el reemplazo de las muestras durante la extracción (bagging) o no (pasting). También se puede tomar una parte de las características en cada subconjunto (Random Patches), o incluso tomar todos los datos y solo una muestra de las características (Random Subspaces). Para el caso de las imágenes del conjunto de partida parece más adecuada la técnica del pasting, no considerando el reemplazo debido a que ya hay imágenes repetidas de algunos diagnósticos. Alternativamente se puede recurrir a un ensemble secuencial (boosting), en el que los estimadores se construyen secuencialmente, intentando reducir el sesgo del estimador final, buscando que la combinación de modelos débiles proporcione un resultado fuerte. En él, cada modelo usa todos los datos de entrenamiento, pero ponderados (remuestreo). Cada iteración se incrementa el peso de los datos mal clasificados para que el siguiente modelo los encuentre más importantes y sea más probable que los clasifique bien. El resultado es una mediana ponderada. En general, en problemas de clasificación y regresión el Ensemble secuencial (boosting) es mejor que en Ensemble de promediación (bagging), pero es más sensible al ruido (presencia de outliers), y tiende al sobreajuste. Sin embargo, no se ha recurrido a ninguna de estas técnicas más complejas porque todas ellas multiplican exponencialmente los tiempos de proceso, y se ha comprobado que con un simple submuestreo los resultados obtenidos no indican que fuese posible obtener mucha mayor mejora con la mayoría de los modelos probados, respecto del conjunto original. A pesar de todo, se usará el resampling simple en la solución final, también por el ahorro de costes que supone reducir el conjunto de datos. Los códigos completos pueden descargarse el repositorio https://github.com/watershed-lab/pneumo-detector En la próxima entrega revisaremos algunos modelos tradicionales para el diagnóstico sobre imágenes.