{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.14","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":71549,"databundleVersionId":8561470,"sourceType":"competition"},{"sourceId":181280574,"sourceType":"kernelVersion"},{"sourceId":200419173,"sourceType":"kernelVersion"}],"dockerImageVersionId":30786,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"Change the file path to match your own environment  ","metadata":{}},{"cell_type":"code","source":"!pip install \"python-gdcm\" pydicom dicomsdl","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-10-12T07:15:34.586018Z","iopub.execute_input":"2024-10-12T07:15:34.587570Z","iopub.status.idle":"2024-10-12T07:15:46.866965Z","shell.execute_reply.started":"2024-10-12T07:15:34.587506Z","shell.execute_reply":"2024-10-12T07:15:46.865232Z"}},"outputs":[],"execution_count":null},{"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)\nimport os\nfrom glob import glob\nimport albumentations as A\nimport torch\nfrom torch.utils.data import Dataset\nfrom torch.utils.data import DataLoader\nimport copy\nfrom scipy.special import softmax\ndef sigmoid(x):\n    return 1/(1 + np.exp(-x))\nfrom sklearn.metrics import roc_auc_score, confusion_matrix, mean_squared_error, average_precision_score, recall_score\nimport warnings\nwarnings.simplefilter('ignore')\n\nimport os\nimport cv2\nimport gdcm\nimport pydicom\nimport zipfile\nimport dicomsdl\nfrom tqdm import tqdm\n\ndef standardize_pixel_array(dcm: pydicom.dataset.FileDataset) -> np.ndarray:\n    \"\"\"\n    Source : https://www.kaggle.com/competitions/rsna-2023-abdominal-trauma-detection/discussion/427217\n    \"\"\"\n    # Correct DICOM pixel_array if PixelRepresentation == 1.\n    pixel_array = dcm.pixel_array\n    if dcm.PixelRepresentation == 1:\n        bit_shift = dcm.BitsAllocated - dcm.BitsStored\n        dtype = pixel_array.dtype \n        new_array = (pixel_array << bit_shift).astype(dtype) >>  bit_shift\n        pixel_array = pydicom.pixel_data_handlers.util.apply_modality_lut(new_array, dcm)\n    return pixel_array\n\ndef get_center_x_path(study_id, series_id):\n    path_x_map = {}\n    for dcm_path in sorted(glob(f\"{dicom_dir}/{study_id}/{series_id}/*.dcm\")):\n        filename = dcm_path.split('/')[-1].replace('.dcm', '')\n    #         if int(instance_number) % 2 == 1:\n    #             continue\n        dicom = dicomsdl.open(dcm_path)\n        pos_x = dicom['ImagePositionPatient'][0]\n        path_x_map[pos_x] = dcm_path\n        xs = []\n        for i, k in enumerate(sorted(path_x_map.keys())):\n            xs.append(k)\n    return path_x_map[xs[len(xs)//2-1]], path_x_map[xs[len(xs)//2]], path_x_map[xs[len(xs)//2+1]]\n\ndef read_sagittal_x_center_dicom(args, verbose=False):\n    study_id, series_id = args\n    ch_imgs = []\n    dcm_path_3 = get_center_x_path(study_id, series_id)\n    for dcm_path in dcm_path_3:\n        dicom = dicomsdl.open(dcm_path)\n        img = dicom.pixelData(storedvalue = True)\n\n        if dicom['PixelRepresentation'] == 1:\n            bit_shift = dicom['BitsAllocated'] - dicom['BitsStored']\n            dtype = img.dtype\n            img = (img << bit_shift).astype(dtype) >>  bit_shift\n        img = img.astype(np.float32)\n\n        intercept = dicom['RescaleIntercept']\n        slope = dicom['RescaleSlope']\n        if (slope is not None) & (intercept is not None):\n            img = img * slope + intercept\n\n        img = (img - img.min()) / (img.max() - img.min() + 1e-6)\n\n        if dicom.PhotometricInterpretation == \"MONOCHROME1\":\n            img = 1 - img\n        img = (img*255.0).astype('uint8')\n\n        ch_imgs.append(img)\n    if verbose:\n        plt.imshow(ch_imgs[1], 'gray')\n        plt.show()\n    img = np.array(ch_imgs).transpose((1,2,0))\n    cv2.imwrite(f'sagittal_all_images/{study_id}___{series_id}.png', img)\n    \ndef select_path_list(lst, N, offset=0, skip=1):\n    if not lst:\n        return [''] * (2 * N + 1)\n    \n    center = (len(lst) // 2) + offset\n    result = [''] * (2 * N + 1)\n    \n    if 0 <= center < len(lst):\n        result[N] = lst[center]\n    \n    for i in range(1, N + 1):\n        left_index = center - i*skip\n        right_index = center + i*skip\n        \n        if 0 <= left_index < len(lst):\n            result[N - i] = lst[left_index]\n        else:\n            result[N - i] = lst[0]\n        if 0 <= right_index < len(lst):\n            result[N + i] = lst[right_index]\n        else:\n            result[N + i] = lst[len(lst)-1]\n    \n    # st()\n    return result\n\n    \ndef read_sagittal_dicom(args, verbose=False):\n    study_id, series_id = args\n    imgs = {}\n    origin_paths = []\n    zs = []\n    xyzs = []\n    paths = []\n    for dcm_path in glob(f\"{dicom_dir}/{study_id}/{series_id}/*.dcm\"):\n        filename = dcm_path.split('/')[-1].replace('.dcm', '')\n#         if int(instance_number) % 2 == 1:\n#             continue\n        dicom = dicomsdl.open(dcm_path)\n        img = dicom.pixelData(storedvalue = True)\n\n        if dicom['PixelRepresentation'] == 1:\n            bit_shift = dicom['BitsAllocated'] - dicom['BitsStored']\n            dtype = img.dtype\n            img = (img << bit_shift).astype(dtype) >>  bit_shift\n        img = img.astype(np.float32)\n\n        intercept = dicom['RescaleIntercept']\n        slope = dicom['RescaleSlope']\n        if (slope is not None) & (intercept is not None):\n            img = img * slope + intercept\n\n        img = (img - img.min()) / (img.max() - img.min() + 1e-6)\n\n        if dicom.PhotometricInterpretation == \"MONOCHROME1\":\n            img = 1 - img\n        img = (img*255.0).astype('uint8')\n        save_path = f'sagittal_all_images/{study_id}___{series_id}___{filename}.png'\n        paths.append(save_path)\n\n        imgs[save_path] = img\n        origin_paths.append(dcm_path)\n        xyzs.append(dicom['ImagePositionPatient'])\n\n    df = pd.DataFrame({\n        'path': paths,\n        'origin_path': origin_paths,\n    })\n    df[['x_pos', 'y_pos', 'z_pos']] = np.array(xyzs)\n    df['study_id'] = study_id\n    df['series_id'] = series_id\n    df['instance_number'] = df['path'].apply(lambda x: int(x.split('___')[-1].replace('.png', '')))\n    df = df.sort_values(['x_pos', 'instance_number'])\n    df = df.drop_duplicates('x_pos')        \n    path_list = df.path.values\n    for path_n, path in enumerate(df.path):\n        if path_n == 0:\n            prev_path = df.path.values[0]\n        else:\n            prev_path = df.path.values[path_n-1]\n        prev_im = imgs[prev_path]            \n        \n        im = imgs[path]\n        if path_n == len(df)-1:\n            next_path = df.path.values[-1]\n        else:\n            next_path = df.path.values[path_n+1]\n        next_im = imgs[next_path]\n  \n        if not (prev_im.shape == im.shape == next_im.shape):\n\n#             print(prev_im.shape, im.shape, next_im.shape)\n            s = prev_im.shape\n            im = cv2.resize(im, s)\n            next_im = cv2.resize(next_im, s)\n\n        image = np.array([prev_im, im, next_im]).transpose((1,2,0))\n        cv2.imwrite(path, image)\n        # if verbose:\n        #     if path_n == 0:\n        #         print(path)\n\n    paths = select_path_list(path_list.tolist(), N=6, offset=0)\n    assert len(paths)==13\n    paths = ','.join(np.array(paths)[[1,2,3,9,10,11]].tolist())\n    df['6ims_paths_123_91011'] = paths\n\n    paths = select_path_list(path_list.tolist(), N=2, offset=0, skip=2)\n    assert len(paths)==5\n    paths = ','.join(paths)\n    df['center_5ims_paths'] = paths\n\n    \n    if verbose:\n        rows = 7\n        for n, k in enumerate(sorted(imgs.keys())):\n            if n % rows == 0:\n                fig = plt.figure(figsize=(20, 20))\n            fig.add_subplot(1, rows, n%rows+1)\n\n            im = imgs[k]\n            plt.imshow(im, 'gray')\n\n            if n % rows == rows-1:\n                plt.show()        \n        \n    return df\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-10-12T07:15:46.869303Z","iopub.execute_input":"2024-10-12T07:15:46.869716Z","iopub.status.idle":"2024-10-12T07:15:51.064044Z","shell.execute_reply.started":"2024-10-12T07:15:46.869673Z","shell.execute_reply":"2024-10-12T07:15:51.062890Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"demo = True","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-10-12T07:17:43.953572Z","iopub.execute_input":"2024-10-12T07:17:43.959228Z","iopub.status.idle":"2024-10-12T07:17:43.971155Z","shell.execute_reply.started":"2024-10-12T07:17:43.959124Z","shell.execute_reply":"2024-10-12T07:17:43.969533Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"df = pd.read_csv('/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_series_descriptions.csv')\nif demo:\n    df = df[df.study_id.isin(df.study_id.unique()[:10])] # Due to lack of storage.\n\naxial_direction = pd.read_csv('/kaggle/input/axial-level-estimation/axial_direction.csv')\ndicom_dir = f'/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_images'","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-10-12T07:17:48.191688Z","iopub.execute_input":"2024-10-12T07:17:48.192119Z","iopub.status.idle":"2024-10-12T07:17:48.778083Z","shell.execute_reply.started":"2024-10-12T07:17:48.192079Z","shell.execute_reply":"2024-10-12T07:17:48.775747Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"image_save_dir = 'axial_all_images'\nos.makedirs(image_save_dir, exist_ok=True)\n\naxial_df = df[df.series_description == 'Axial T2']\nargs = axial_df.drop_duplicates(['study_id', 'series_id'])[['study_id', 'series_id']].values\n\n\ndef read_axial_dicom(args, verbose=False):\n    study_id, series_id = args\n    imgs = {}\n    origin_paths = []\n    zs = []\n    xyzs = []\n    paths = []\n    series_axial_direction = axial_direction[axial_direction.series_id==series_id]\n    for dcm_path in glob(f\"{dicom_dir}/{study_id}/{series_id}/*.dcm\"):\n        filename = dcm_path.split('/')[-1].replace('.dcm', '')\n        dicom = dicomsdl.open(dcm_path)\n        img = dicom.pixelData(storedvalue = True)\n\n        if dicom['PixelRepresentation'] == 1:\n            bit_shift = dicom['BitsAllocated'] - dicom['BitsStored']\n            dtype = img.dtype\n            img = (img << bit_shift).astype(dtype) >>  bit_shift\n        img = img.astype(np.float32)\n\n        intercept = dicom['RescaleIntercept']\n        slope = dicom['RescaleSlope']\n        if (slope is not None) & (intercept is not None):\n            img = img * slope + intercept\n\n        pos_z = dicom['ImagePositionPatient'][-1]\n        \n        img = (img - img.min()) / (img.max() - img.min() + 1e-6)\n\n        if dicom.PhotometricInterpretation == \"MONOCHROME1\":\n            img = 1 - img\n        img = (img*255.0).astype('uint8')\n        save_path = f'axial_all_images/{study_id}___{series_id}___{filename}.png'\n        paths.append(save_path)\n\n        imgs[save_path] = img\n        origin_paths.append(dcm_path)\n        xyzs.append(dicom['ImagePositionPatient'])\n\n    df = pd.DataFrame({\n        'path': paths,\n        'origin_path': origin_paths,\n    })\n    df[['x_pos', 'y_pos', 'z_pos']] = np.array(xyzs)\n    df['study_id'] = study_id\n    df['series_id'] = series_id\n    df['instance_number'] = df['path'].apply(lambda x: int(x.split('___')[-1].replace('.png', '')))\n\n    l = len(df)\n    if len(df.merge(series_axial_direction[['instance_number', 'z']], on='instance_number')) == l:\n        df = df.merge(series_axial_direction[['instance_number', 'z']], on='instance_number')\n        df = df.sort_values(['z', 'instance_number'])\n    else:\n        df = df.sort_values(['z_pos', 'instance_number'])\n    path_list = df.path.values\n\n    for path_n, path in enumerate(df.path):\n        im = imgs[path]\n        if path_n == 0:\n            prev_path = df.path.values[0]\n        else:\n            prev_path = df.path.values[path_n-1]\n        prev_im = imgs[prev_path]\n        \n        if path_n == len(df)-1:\n            next_path = df.path.values[-1]\n        else:\n            next_path = df.path.values[path_n+1]\n        next_im = imgs[next_path]\n  \n        if not (prev_im.shape == im.shape == next_im.shape):\n\n            s = prev_im.shape\n            im = cv2.resize(im, s)\n            next_im = cv2.resize(next_im, s)\n\n        image = np.array([prev_im, im, next_im]).transpose((1,2,0))\n        cv2.imwrite(path, image)\n        \n    if verbose:\n        rows = 7\n        for n, k in enumerate(sorted(imgs.keys())):\n            if n % rows == 0:\n                fig = plt.figure(figsize=(20, 20))\n            fig.add_subplot(1, rows, n%rows+1)\n\n            im = imgs[k]\n            plt.imshow(im, 'gray')\n\n            if n % rows == rows-1:\n                plt.show()        \n        \n    return df\n\nfrom multiprocessing import Pool\n\np = Pool(processes=4)\nresults = []\nwith tqdm(total=len(args)) as pbar:\n    for res in p.imap(read_axial_dicom, args):\n        results.append(res)\n        pbar.update(1)\np.close()\n\naxial_df = pd.concat(results)\naxial_df.to_csv('axial_df.csv', index=False)\n\nimage_save_dir = 'sagittal_all_images'\nos.makedirs(image_save_dir, exist_ok=True)\n\nsagittal_df = df[df.series_description != 'Axial T2']\nargs = sagittal_df.drop_duplicates(['study_id', 'series_id'])[['study_id', 'series_id']].values\n\nfrom multiprocessing import Pool\n\np = Pool(processes=4)\nresults = []\nwith tqdm(total=len(args)) as pbar:\n    for res in p.imap(read_sagittal_dicom, args):\n        results.append(res)\n        pbar.update(1)\np.close()\n\nsagittal_df = pd.concat(results)\nsagittal_df.to_csv('sagittal_df.csv', index=False)\n\nsagittal_df = pd.read_csv('sagittal_df.csv')\nseries_description_df = pd.read_csv(f'/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_series_descriptions.csv')\nsagittal_df = sagittal_df.merge(series_description_df, on=['study_id', 'series_id'])\n\nfor study_id_n, (study_id, idf) in enumerate(tqdm(sagittal_df.groupby('study_id'))):\n    t1 = idf[idf.series_description=='Sagittal T1'].sort_values('instance_number')\n    if len(t1)!=0:\n        t1 = t1[t1.series_id==t1.series_id.values[0]]\n    t2 = idf[idf.series_description=='Sagittal T2/STIR'].sort_values('instance_number')\n    if len(t2)!=0:\n        t2 = t2[t2.series_id==t2.series_id.values[0]]\n    if len(t1)==0:\n        m = t2.instance_number.max()\n        mi = t2.instance_number.min()\n    elif (len(t2)==0):\n        m = t1.instance_number.max()\n        mi = t1.instance_number.min()\n    else:\n        m = max([t1.instance_number.max(), t2.instance_number.max()])\n        mi = min([t1.instance_number.min(), t2.instance_number.min()])\n    for n in range(mi, m+1):\n        n1 = t1[t1.instance_number == n]\n        n2 = t2[t2.instance_number == n]\n        if len(n1) != 0:\n            # print(n1.path.values[0])\n            im1 = cv2.imread(n1.path.values[0])[:,:,1]\n        else:\n            im1 = None\n        if len(n2) != 0:\n            # print(n2.path.values[0])\n            im2 = cv2.imread(n2.path.values[0])[:,:,1]\n        else:\n            im2 = None\n        if im1 is None:\n            # raise\n            im1 = np.zeros(im2.shape)\n        if im2 is None:\n            # raise\n            im2 = np.zeros(im1.shape)\n        if im1.shape!=im2.shape:\n            im1 = cv2.resize(im1, (im2.shape[1], im2.shape[0]))\n        im = np.array([im1, im2, im1]).transpose((1,2,0))\n\n        cv2.imwrite(f'sagittal_all_images/{study_id}___{n}.png', im)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-10-12T07:17:48.781103Z","iopub.execute_input":"2024-10-12T07:17:48.782004Z","iopub.status.idle":"2024-10-12T07:18:08.985172Z","shell.execute_reply.started":"2024-10-12T07:17:48.781942Z","shell.execute_reply":"2024-10-12T07:18:08.983581Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"axial_df[['study_id', 'series_id', 'path', 'instance_number']]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-10-12T07:18:08.988210Z","iopub.execute_input":"2024-10-12T07:18:08.989303Z","iopub.status.idle":"2024-10-12T07:18:09.024450Z","shell.execute_reply.started":"2024-10-12T07:18:08.989239Z","shell.execute_reply":"2024-10-12T07:18:09.022707Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}