DETECCIÓN DE NEUMONÍA RADIOGRAFÍAS: 1 – CARGA Y EXPLORACIÓN DE DATOS

Publicado por @JoseMa_deCuenca el March 8, 2020, 9:10 p.m.
El uso de la visión por computadora para diagnóstico médico es un clásico entre los problemas de inteligencia artificial. Aunque hasta hace poco era muy difícil abordar desde un ordenador personal. Con el incremento de la potencia de los procesadores, la aparición de nuevas librerías, y sobre todo, teniendo en cuenta algunas buenas prácticas durante la programación, hoy es posible obtener elevados rendimientos y resolver problemas complejos con modelos que pueden entrenarse en solo unas horas con nuestro propio PC. Vamos a ver como se aplicaría a un caso real, a través de varias entregas: el diagnóstico de pneumonía a partir de radiografías pulmonares. El ejemplo puede parecer escogido por el despliegue mediático y la histeria colectiva que está desatando la crisis del coronavirus, pero no es así. Fue desarrollado medio año antes del descubrimiento del virus, y elegido como trabajo de fin de máster porque pertenecía a un campo que me resulta totalmente ajeno, percísamente para probar la bondad de los conocimientos aprendidos sin ninguna sinergia o cruce con experiencia previa. En general, un problema de este tipo se puede tratar como un clasificador de imágenes utilizando aprendizaje supervisado. Para resolverlo con éxito, los datos de partida deben ser: - Correctos, con un etiquetado perfecto entre las diversas categorías, generalmente enfermos y sanos, aunque también se pueden discriminar causas, gravedad, tipos de enfermedades… - Numerosos, para favorecer la generalización e incluso permitir selección de subconjuntos aleatorios. - Datos generales, con variedad suficiente para evitar sesgos que puedan pasar al modelo durante el aprendizaje. - Relevantes, con las características necesarias y solo con ellas para evitar complejidades accesorias. - Suficientes con información representativa para recoger todos los patrones que debe detectar la herramienta de aprendizaje. - Actuales, de manera que recojan las variaciones temporales en caso de producirse, especialmente online. - Equilibrados entre las diferentes categorías, para que el aprendizaje de todos los tipos sea homogéneo. En realidad el dataset de origen está sesgado, pero se corrige más adelante. Para entrenar el modelo diagnosticador de pneumonía a partir de radiografías usaremos un conjunto de imágenes libres, el dataset disponible en kaggle con más de 5.000 radiografías de tórax, clasificadas en 2 categorías según pertenezcan a enfermos de pneumonía o no (https://www.kaggle.com/paultimothymooney/chest_xray_pneumonia), aunque debido al interés que puede tener identificar el origen de la causa de la neumonía, generaremos también una clasificación que discrimine si es bacteriana o vírica. Llegar a saber ya si fue causada por el COVID-19 con estos modelos hoy no es posible, pero tampoco me atreveo a decir que resulte imposible: quien sabe si los artículos llegan a inspirar a alguien. Vamos con ello. Los pasos a seguir los recorreremos siguiendo el esquema, en varios artículos. ![enter image description here][1] La forma más interesante de procesar los datos es a partir de archivos individuales alojados en disco. Cualquier intento de pasarles a matriz numérica conlleva unos recursos que imposibilitarían el uso de un ordenador personal, o al menos se ralentizaría mucho todo el proceso. Esto se produce básicamente porque el tiempo de carga de un gran archivo es muy elevado frente al de lectura de varios miles de ficheros mucho más livianos. Se puede comprobar con el siguiente código, con el que además es un generador de tabla de características, que nos facilitará enormemente el análisis exploratorio de los datos. Incluso genera una clasificación adicional a la del dataset de partida, en 3 clases, mediante el procesado de los nombres de archivo. # DATAFRAME DE ARCHIVOS DE ENTRADA ORIGINALES %matplotlib inline import os import numpy as np import pandas as pd from PIL import Image # Uso el sistema para obtener la Lista con todos los ficheros de la ruta path = './data/chest_xray' # Variable para la ruta al directorio lstDir = os.walk(path) #os.walk()Lista directorios y ficheros # Creo lista para extraer los datos de la anterior del sistema # Lista vacia para rellenar con datos de los ficheros lstFiles = [] for root, dirs, files in lstDir: for fichero in files: # Separo nombre del fichero de la extension (nombreFichero, extension) = os.path.splitext(fichero) # Controlo que solo se procesen los ficheros jpeg (evitando los .DS_Store u otros) if extension == ".jpeg": # Compongo la ruta completa, convirtiendo a tipo Unix para independizar de SO. Ojo escape ruta = root.replace("\\", "/")+"/"+nombreFichero+extension # Uso la ruta completa como identificador de la imagen lstFiles.append(ruta) # Extraigo el diagnostico if root.find("NORMAL")>0: lstFiles.append("NORMAL") else: lstFiles.append("PNEUMONIA") # Extraigo el tipo de infeccion if root.find("PNEUMONIA")>0: if nombreFichero.find("bacteria")>0: lstFiles.append("BACTERIA") else: lstFiles.append("VIRUS") else: lstFiles.append("NONE") # Cargo cada imagen para sacar caracteristicas internas imagen = Image.open(ruta) lstFiles.append(imagen.mode) lstFiles.append(imagen.size[0]) # anchura lstFiles.append(imagen.size[1]) # altura lstFiles.append(round(imagen.size[0]*imagen.size[1]/1000000, 2)) # resolucion Megapixeles lstFiles.append(round(imagen.size[0]/imagen.size[1], 4)) # relacion de aspecto # Convierto la segunda lista en array del sistema npdf_files = np.array(lstFiles) # Lo llevo al formato de columnas ordenadas npdf_files = npdf_files.reshape([-1,8]) # Compruebo que todo es correcto print(npdf_files.shape) print(npdf_files) # Creo un dataframe de pandas df_file = pd.DataFrame(npdf_files) # Lo grabo como csv df_file.to_csv (path+'/df_file.csv', index = None, header=False) # Compruebo visualmente df_file.head(10) Descartar el uso de una matriz numérica como dataset no es obstáculo para comenzar realizando un análisis exploratorio lo más exhaustivo posible de los datos de partida. Para ello podemos utilizar la tabla de características que acabamos de generar: ## ANALISIS ESTADISTICO import matplotlib.pyplot as plt # Tabla cruzada de casos tab = pd.crosstab(df_file[1], df_file[2]) print(tab) # Grafico de registros por categoria # Columna de diagnostico df_file[1].value_counts().plot(kind='bar', alpha=0.5) plt.title("Diagnóstico") plt.show() # Columna de tipo infección df_file[2].value_counts().plot(kind='bar', alpha=0.5) plt.title("Causa") plt.show() # Tipo de imagen (codigos de Pil image.mode) # 1 (1-bit pixels, black and white, stored with one pixel per byte) # L (8-bit pixels, black and white) # P (8-bit pixels, mapped to any other mode using a colour palette) # RGB (3x8-bit pixels, true colour) # RGBA (4x8-bit pixels, true colour with transparency mask) # CMYK (4x8-bit pixels, colour separation) # YCbCr (3x8-bit pixels, colour video format) # I (32-bit signed integer pixels) # F (32-bit floating point pixels) # Columna de tipo imagen df_file[3].value_counts().plot(kind='bar', alpha=0.5) plt.title("Tipo de imagen") plt.show() # Indago las imagenes que estan en color por el tipo infeccion tab = pd.crosstab(df_file[2], df_file[3]) print(tab) ## Análisis de imágenes por rango de resolucion # Defino rangos de resolucion en MegaPixeles rangos_resol = ['0-1', '1-2', '2-3', '3-4', '>4'] # Cuento el numero de registros en cada rango reso1 = df_file[(df_file[6].astype(float)< 1.0)][6].count() reso2 = df_file[(df_file[6].astype(float)< 2.0)][6].count() - reso1 reso3 = df_file[(df_file[6].astype(float)< 3.0)][6].count() - reso2 - reso1 reso4 = df_file[(df_file[6].astype(float)< 4.0)][6].count() - reso3 - reso2 - reso1 reso5 = df_file[(df_file[6].astype(float)>= 4.0)][6].count() # Creo lista de valores por rangos resol = [reso1, reso2, reso3, reso4, reso5] # Represento con gráfico de barras plt.bar(rangos_resol, resol) plt.xlabel("Resolucion") plt.ylabel("Nº imágenes") plt.title("Imágenes por rango de resolucion") plt.show() ## Análisis de imágenes por relación de aspecto # Defino rangos de resolucion en MegaPixeles rangos_aspect = ['<1.0', '1-1.25', '1.25-1.35', '1.35-1.5', '1.5-1.65' , '>1.65'] # Cuento el número de registros en cada rango aspe1 = df_file[(df_file[7].astype(float)< 1.0)][7].count() aspe2 = df_file[(df_file[7].astype(float)< 1.25)][7].count() - aspe1 aspe3 = df_file[(df_file[7].astype(float)< 1.35)][7].count() - aspe2 - aspe1 aspe4 = df_file[(df_file[7].astype(float)< 1.5)][7].count() - aspe3 - aspe2 - aspe1 aspe5 = df_file[(df_file[7].astype(float)< 1.65)][7].count() - aspe4 - aspe3 - aspe2 - aspe1 aspe6 = df_file[(df_file[7].astype(float)>= 1.65)][7].count() # Creo lista de valores por rangos aspect = [aspe1, aspe2, aspe3, aspe4, aspe5, aspe6] # Represento con gráfico de barras plt.bar(rangos_aspect, aspect) plt.xlabel("Relación de aspecto") plt.ylabel("Nº imágenes") plt.title("Imágenes por relación de aspecto") plt.show() A pesar de descartar el uso de un único gran archivo de matriz numérica durante el preprocesado de las imágenes, dado su elevado número, debemos ser cuidadosos a la hora de seleccionar el método más apropiado para su carga. Probaremos dos métodos con el conjunto de validación, que es el más reducido (16 elementos en esa carpeta), para demostrar el tiempo de proceso que podemos ahorrar: básicamente consisten en su redimensionamiento y conversión a np. # PRUEBA PARA MEJORAR TIEMPOS DE CARGA DE DATOS # Preparo conjuntos de datos # Creo array con imagenes para modelo # Creo listas con las rutas de las imágenes para cada conjunto VAL_files = (df_file[df_file['conjunto']=='VAL']['fichero'].iloc[:]).tolist() # Creo array con etiquetas para modelo, por ejemplo con columna diagnostico categorizada. Si patogeno, cambiar VAL_labels = (df_file[df_file['conjunto']=='VAL']['diag_ct'].iloc[:]).tolist() # Creo listas para cargar los datos de las imagenes y las etiquetas de cada conjunto X_val = [] y_val = [] # Comienzo cargando los datos del conjunto de validacion por ser mas reducido # Pruebo primero con resize de sckit-learn print("Procesando dataset de datos:", len(VAL_labels),"elementos") # Inicia cronometro start_time = time.time() for img in VAL_files: img = cv2.imread(str(img)) img = resize(img, (150, 150, 3)) img = np.asarray(img) X_val.append(img) X_val = np.array(X_val) # Continua cargando los datos de las etiquetas for label in VAL_labels: y_val.append(label) y_val = np.array(y_val) # Detiene el cronometro end_time = time.time() # Informa cronometro dt = round(end_time - start_time, 2) print("Este conjunto tardó en cargarse: ", dt, "segundos (", round(dt/60,1) , "minutos)") # Informa de la carga print(X_val.shape) print(y_val.shape) # PRUEBA ALTERNATIVA PARA MEJORAR TIEMPO DE CARGA DE DATOS # Creo listas con las rutas de las imágenes para cada conjunto VAL_files = (df_file[df_file['conjunto']=='VAL']['fichero'].iloc[:]).tolist() # Creo array con etiquetas para modelo, por ejemplo con columna diagnostico categorizada. Si patogeno, cambiar VAL_labels = (df_file[df_file['conjunto']=='VAL']['diag_ct'].iloc[:]).tolist() # Creo listas para cargar los datos de las imagenes y las etiquetas de cada conjunto X_val = [] y_val = [] # Comienzo cargando los datos del conjunto de validacion por ser mas reducido print("Procesando dataset de datos:", len(VAL_labels),"elementos") # Inicia cronometro start_time = time.time() for img in VAL_files: img = cv2.imread(str(img)) # img = resize(img, (150, 150, 3)) # Alternativamente img = cv2.resize(img, (150, 150)) if img.shape[2] ==1: img = np.dstack([img, img, img]) img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB) # Cambio del formato de entrada (orden BGR de OpenCV) al estandar (orden RGB) # Se puede probar con otros espacios de color # img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) # Cambia de entrada (BGR OpenCV) a escala de grises # img = cv2.cvtColor(img, cv2.COLOR_BGR2HSV) # Cambia de entrada (BGR OpenCV) a HSV (matiz, saturación, iluminancia) img = img.astype(np.float32)/255.0 # aprovecho para normalizar los datos. Se corresponde a la escala del espacio de color img = np.asarray(img) X_val.append(img) X_val = np.array(X_val) # Continua cargando los datos de las etiquetas for label in VAL_labels: y_val.append(label) y_val = np.array(y_val) # Detiene el cronometro end_time = time.time() # Informa cronometro dt = round(end_time - start_time, 2) print("Este conjunto tardó en cargarse: ", dt, "segundos (", round(dt/60,1) , "minutos)") # Informa de la carga print(X_val.shape) print(y_val.shape) Podemos comprobar que es mucho más eficiente el segundo método, que realiza un redimensionamiento solo en resolución, no en canales; realiza una transformación a RGB estándar (CV2 trabaja sobre BGR) y normaliza los datos antes de su conversión a una matriz numpy. En un i7 8700 con 16Gb de RAM este método logra reducir la carga de las 16 imágenes de 2,17sg a solo 0,17sg (doce veces más rápido). Hasta aquí la primera parte, seguiremos en breve para ver cómo se usa la transformación y segmentación de las imágenes. Los códigos completos pueden descargarse el repositorio https://github.com/watershed-lab/pneumo-detector [1]: https://raw.githubusercontent.com/watershed-lab/pneumo-detector/master/esquema.png