{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"seed = 2022","metadata":{"execution":{"iopub.status.busy":"2022-11-21T10:52:45.12234Z","iopub.execute_input":"2022-11-21T10:52:45.122733Z","iopub.status.idle":"2022-11-21T10:52:45.128394Z","shell.execute_reply.started":"2022-11-21T10:52:45.122702Z","shell.execute_reply":"2022-11-21T10:52:45.127102Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\nfrom sklearn.model_selection import train_test_split\n\nimport os\n\nimport seaborn as sns\nimport matplotlib\nimport matplotlib.pyplot as plt\n\nimport pydicom\nimport pylab","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-11-21T10:52:46.061914Z","iopub.execute_input":"2022-11-21T10:52:46.062329Z","iopub.status.idle":"2022-11-21T10:52:46.068953Z","shell.execute_reply.started":"2022-11-21T10:52:46.062299Z","shell.execute_reply":"2022-11-21T10:52:46.067218Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train = pd.read_csv(\"../input/rsna-intracranial-hemorrhage-detection/rsna-intracranial-hemorrhage-detection/stage_2_train.csv\")\n\ntrain[\"type\"] = train[\"ID\"].str.split(\"_\", expand = True)[2]\ntrain[\"ID\"] = \"ID_\" + train[\"ID\"].str.split(\"_\", expand = True)[1]\n\ntrain = train.pivot_table(index = \"ID\", columns = \"type\", values = \"Label\")\n\n# Quitar imagen corrupta de la BDD\ntrain = train.drop([\"ID_6431af929\"], axis = 0)","metadata":{"execution":{"iopub.status.busy":"2022-11-21T10:56:29.796731Z","iopub.execute_input":"2022-11-21T10:56:29.797869Z","iopub.status.idle":"2022-11-21T10:57:12.038769Z","shell.execute_reply.started":"2022-11-21T10:56:29.797826Z","shell.execute_reply":"2022-11-21T10:57:12.037524Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_images_path = \"../input/rsna-intracranial-hemorrhage-detection/rsna-intracranial-hemorrhage-detection/stage_2_train/\"","metadata":{"execution":{"iopub.status.busy":"2022-11-21T10:57:12.0406Z","iopub.execute_input":"2022-11-21T10:57:12.041083Z","iopub.status.idle":"2022-11-21T10:57:12.046173Z","shell.execute_reply.started":"2022-11-21T10:57:12.041047Z","shell.execute_reply":"2022-11-21T10:57:12.044951Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mkdir train_images","metadata":{"execution":{"iopub.status.busy":"2022-11-21T10:53:36.565985Z","iopub.execute_input":"2022-11-21T10:53:36.567293Z","iopub.status.idle":"2022-11-21T10:53:37.668682Z","shell.execute_reply.started":"2022-11-21T10:53:36.567231Z","shell.execute_reply":"2022-11-21T10:53:37.667363Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ls","metadata":{"execution":{"iopub.status.busy":"2022-11-21T10:54:02.301183Z","iopub.execute_input":"2022-11-21T10:54:02.301602Z","iopub.status.idle":"2022-11-21T10:54:03.405652Z","shell.execute_reply.started":"2022-11-21T10:54:02.301566Z","shell.execute_reply":"2022-11-21T10:54:03.404352Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_first_of_dicom_field_as_int(x):\n    #get x[0] as in int is x is a 'pydicom.multival.MultiValue', otherwise get int(x)\n    if type(x) == pydicom.multival.MultiValue:\n        return int(x[0])\n    else:\n        return int(x)\n\ndef get_windowing(data):\n    dicom_fields = [data[('0028','1050')].value, #window center\n                    data[('0028','1051')].value, #window width\n                    data[('0028','1052')].value, #intercept\n                    data[('0028','1053')].value] #slope\n    return [get_first_of_dicom_field_as_int(x) for x in dicom_fields]\n\ndef window_image(img, window_center,window_width, intercept, slope):\n    img = (img*slope +intercept)\n    img_min = window_center - window_width//2\n    img_max = window_center + window_width//2\n    img[img<img_min] = img_min\n    img[img>img_max] = img_max\n    return img ","metadata":{"execution":{"iopub.status.busy":"2022-11-21T10:57:12.048084Z","iopub.execute_input":"2022-11-21T10:57:12.048528Z","iopub.status.idle":"2022-11-21T10:57:12.061937Z","shell.execute_reply.started":"2022-11-21T10:57:12.048485Z","shell.execute_reply":"2022-11-21T10:57:12.060609Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_borders(img):\n    x1, y1 = 0, 0\n    y2, x2 = img.shape\n\n    for i in range(0, len(img[0])-1,1):\n        if np.sum(img[:,i]) != np.sum(img[:,i+1]):\n            x1 = i\n            break\n\n    for i in range(len(img[0])-2,0,-1):\n        if np.sum(img[:,i]) != np.sum(img[:,i+1]):\n            x2 = i\n            break\n        \n    for j in range(0,len(img)-1, 1):    \n        if np.sum(img[j])!=np.sum(img[j+1]):\n            y1 = j\n            break\n        \n    for j in range(len(img)-2, 0,-1):\n        if np.sum(img[j])!=np.sum(img[j+1]):\n            y2 = j\n            break\n    \n    return x1, x2, y1, y2","metadata":{"execution":{"iopub.status.busy":"2022-11-21T10:57:12.065415Z","iopub.execute_input":"2022-11-21T10:57:12.065939Z","iopub.status.idle":"2022-11-21T10:57:12.07927Z","shell.execute_reply.started":"2022-11-21T10:57:12.065891Z","shell.execute_reply":"2022-11-21T10:57:12.077887Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def normalize(arr):\n    arr_min = np.min(arr)\n    arr_max = np.max(arr)\n    if arr_max-arr_min == 0:\n        return np.zeros(arr.shape)\n    else:\n        return (arr-arr_min)/(arr_max-arr_min)  ","metadata":{"execution":{"iopub.status.busy":"2022-11-21T10:57:17.341137Z","iopub.execute_input":"2022-11-21T10:57:17.341566Z","iopub.status.idle":"2022-11-21T10:57:17.34783Z","shell.execute_reply.started":"2022-11-21T10:57:17.34153Z","shell.execute_reply":"2022-11-21T10:57:17.346896Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from skimage.filters import laplace\nfrom skimage.filters import gaussian\nimport cv2\nfrom skimage.transform import resize\n\n# centra la imagen, quitando el sobrante y ajustando al tamaño original\n# ret es un parametro para ver los pasos realizados a las imagenes\ndef image_adjust(path, name, h, w, ret = False, save_dir = None):\n    ds = pydicom.dcmread(path + name + \".dcm\")\n    \n    # corregir windowing\n    image = ds.pixel_array\n    window_center , window_width, intercept, slope = get_windowing(ds)\n    image_windowed = window_image(image, window_center, window_width, intercept, slope)\n    \n    # hay que hacer un marco exterior, para que si hay imagenes que están cortadas se puede dibujar el contorno exterior\n    image_windowed[0] = np.min(image_windowed)\n    image_windowed[:,0] = np.min(image_windowed)\n    image_windowed[image_windowed.shape[0]-1] = np.min(image_windowed)\n    image_windowed[:,image_windowed.shape[1]-1] = np.min(image_windowed)\n    \n    # binarizar la imagen    \n    binary = np.zeros(image_windowed.shape).astype(np.uint8)\n    binary[image_windowed>np.min(image_windowed)] = 255\n    binary[image_windowed==np.min(image_windowed)] = 0\n\n    # aplicar filtro laplaciano\n    lap = laplace(binary).astype(np.uint8).copy()\n\n    # se busca los contornos de la imagen con el filtro laplaciano\n    contours, hierarchy = cv2.findContours(lap.copy(), cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE)\n    \n    # se seleccionan solo los 10 contornos con mas puntos para optimizar los calculos de las máscaras\n    contours_shape = [c.shape[0] for c in contours]\n    contours_shape_sorted = contours_shape.copy()\n    contours_shape_sorted.sort(reverse=True)\n    max_contours = contours_shape_sorted[:20]\n    best_contours = set([contours_shape.index(c) for c in max_contours])\n    # elegir la mejor mascara\n    final_mask = np.zeros(image_windowed.shape)\n    \n    for i in best_contours:\n        # rellenar el contorno de la imagen con el filtro laplaciano para obtener la mascara\n        mask = cv2.drawContours(lap.copy(), contours, contourIdx = i, hierarchy = hierarchy, color=(255, 255, 255), thickness=cv2.FILLED)\n        mask[mask<255] = 0\n        \n        # aplicar filtro gaussiano para borrar pequeñas regiones de la mascara\n        mask = gaussian(mask)\n        mask[mask<1] = 0\n        mask[mask == 1] = 255\n        # comprobar la región dentro de la mascara que más tamaño tiene\n        if np.count_nonzero(mask) > np.count_nonzero(final_mask):\n            final_mask = mask\n\n    # aplicar la mascara resultante a la imagen de partida\n    result = image_windowed.copy()\n    result[final_mask<np.max(final_mask)] = np.min(image_windowed)\n\n    # normalizar la imagen entre 0 y 1\n    normalized = normalize(result)\n    \n    # obtención de bordes para el recorte\n    x1, x2, y1, y2 = get_borders(normalized)\n    \n    # se recorta la imagen con los bordes obtenidos\n    cropped = normalized[y1:y2, x1:x2]\n    \n    # devuelve la imagen con el tamaño deseado\n    resized = resize(cropped, (h, w))\n    \n    # comprobar si se quiere visualizar el proceso\n    if ret:\n        border = normalized.copy()\n        border[y1] = 1\n        border[y2] = 1\n        border[:,x1] = 1\n        border[:,x2] = 1\n        return (image, image_windowed, binary, lap, final_mask, result, border, cropped, resized)\n    \n    if save_dir == None:\n        return resized\n    else:\n        np.save(save_dir + \"/\" + name + \".npy\", resized)","metadata":{"execution":{"iopub.status.busy":"2022-11-21T11:14:24.113306Z","iopub.execute_input":"2022-11-21T11:14:24.113735Z","iopub.status.idle":"2022-11-21T11:14:24.134076Z","shell.execute_reply.started":"2022-11-21T11:14:24.113699Z","shell.execute_reply":"2022-11-21T11:14:24.133056Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_index = train.index.values","metadata":{"execution":{"iopub.status.busy":"2022-11-21T11:14:24.371358Z","iopub.execute_input":"2022-11-21T11:14:24.371969Z","iopub.status.idle":"2022-11-21T11:14:24.376291Z","shell.execute_reply.started":"2022-11-21T11:14:24.371934Z","shell.execute_reply":"2022-11-21T11:14:24.375458Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from joblib import Parallel, delayed","metadata":{"execution":{"iopub.status.busy":"2022-11-21T11:14:24.898991Z","iopub.execute_input":"2022-11-21T11:14:24.899687Z","iopub.status.idle":"2022-11-21T11:14:24.904523Z","shell.execute_reply.started":"2022-11-21T11:14:24.899647Z","shell.execute_reply":"2022-11-21T11:14:24.903447Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"save_path_train = \"./train_images\"","metadata":{"execution":{"iopub.status.busy":"2022-11-21T11:14:25.410307Z","iopub.execute_input":"2022-11-21T11:14:25.410701Z","iopub.status.idle":"2022-11-21T11:14:25.415906Z","shell.execute_reply.started":"2022-11-21T11:14:25.410668Z","shell.execute_reply":"2022-11-21T11:14:25.414465Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nnp.array(Parallel(n_jobs=-1)(delayed(image_adjust)(train_images_path, train.index[i], 256, 256, save_dir = save_path_train) for i in range(100)))","metadata":{"execution":{"iopub.status.busy":"2022-11-21T11:14:26.689545Z","iopub.execute_input":"2022-11-21T11:14:26.690413Z","iopub.status.idle":"2022-11-21T11:14:35.97554Z","shell.execute_reply.started":"2022-11-21T11:14:26.690372Z","shell.execute_reply":"2022-11-21T11:14:35.974109Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"image = np.load(save_path_train + \"/ID_000716c43.npy\")","metadata":{"execution":{"iopub.status.busy":"2022-11-21T11:07:21.234676Z","iopub.execute_input":"2022-11-21T11:07:21.235123Z","iopub.status.idle":"2022-11-21T11:07:21.24252Z","shell.execute_reply.started":"2022-11-21T11:07:21.235087Z","shell.execute_reply":"2022-11-21T11:07:21.241117Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!du -sh","metadata":{"execution":{"iopub.status.busy":"2022-11-21T11:09:53.119304Z","iopub.execute_input":"2022-11-21T11:09:53.119706Z","iopub.status.idle":"2022-11-21T11:09:54.248267Z","shell.execute_reply.started":"2022-11-21T11:09:53.119675Z","shell.execute_reply":"2022-11-21T11:09:54.246657Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!zip -r train_images.zip ./train_images","metadata":{"execution":{"iopub.status.busy":"2022-11-21T11:13:09.774296Z","iopub.execute_input":"2022-11-21T11:13:09.774749Z","iopub.status.idle":"2022-11-21T11:13:12.888935Z","shell.execute_reply.started":"2022-11-21T11:13:09.774714Z","shell.execute_reply":"2022-11-21T11:13:12.887885Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!rm -r train_images","metadata":{"execution":{"iopub.status.busy":"2022-11-21T11:22:40.13213Z","iopub.execute_input":"2022-11-21T11:22:40.132585Z","iopub.status.idle":"2022-11-21T11:22:41.289057Z","shell.execute_reply.started":"2022-11-21T11:22:40.132546Z","shell.execute_reply":"2022-11-21T11:22:41.28719Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}