{"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":"# Installs:\n#!pip install -U pylibjpeg[all]\n#!pip install cython\n#! pip install pylibjpeg pylibjpeg-libjpeg pydicom\n#!pip install -U python-gdcm\n#!pip install pillow\n'''\nimport time\n\n!wget 'https://anaconda.org/conda-forge/libjpeg-turbo/2.1.0/download/linux-64/libjpeg-turbo-2.1.0-h7f98852_0.tar.bz2' -q\n!wget 'https://anaconda.org/conda-forge/libgcc-ng/9.3.0/download/linux-64/libgcc-ng-9.3.0-h2828fa1_19.tar.bz2' -q\n!wget 'https://anaconda.org/conda-forge/gdcm/2.8.9/download/linux-64/gdcm-2.8.9-py37h500ead1_1.tar.bz2' -q\n!wget 'https://anaconda.org/conda-forge/conda/4.10.1/download/linux-64/conda-4.10.1-py37h89c1867_0.tar.bz2' -q\n!wget 'https://anaconda.org/conda-forge/certifi/2020.12.5/download/linux-64/certifi-2020.12.5-py37h89c1867_1.tar.bz2' -q\n!wget 'https://anaconda.org/conda-forge/openssl/1.1.1k/download/linux-64/openssl-1.1.1k-h7f98852_0.tar.bz2' -q\n\n!conda install 'libjpeg-turbo-2.1.0-h7f98852_0.tar.bz2' -c conda-forge -y\n!conda install 'libgcc-ng-9.3.0-h2828fa1_19.tar.bz2' -c conda-forge -y\n!conda install 'gdcm-2.8.9-py37h500ead1_1.tar.bz2' -c conda-forge -y\n!conda install 'conda-4.10.1-py37h89c1867_0.tar.bz2' -c conda-forge -y\n!conda install 'certifi-2020.12.5-py37h89c1867_1.tar.bz2' -c conda-forge -y\n!conda install 'openssl-1.1.1k-h7f98852_0.tar.bz2' -c conda-forge -y\n\n!pip install pylibjpeg pylibjpeg-libjpeg pylibjpeg-openjpeg\n\n# Relevant imports:\nimport scipy\nimport scipy.interpolate as interpolate\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nfrom shutil import copy, make_archive\nimport pydicom\nimport matplotlib.pyplot as plt\nimport cv2\nimport SimpleITK as sitk\nimport os\n\n# We create a function that takes a CT and makes its voxels isotropic:\ndef do_interpolate(image_data, steps):\n    dx, dy, dz = 1., 1., 1.   # Desired (isotropic) spacing\n    image_data = np.moveaxis(image_data, 0, -1)\n    x, y, z = [steps[k] * np.arange(image_data.shape[k]) for k in range(len(steps))]  # original grid\n    # We create an interpolator:\n    f_scan = interpolate.RegularGridInterpolator((x, y, z), image_data, bounds_error = False, method = \"nearest\") \n    # New grid (with adequate spacing):   \n    new_grid = np.mgrid[0:x[-1]:dx, 0:y[-1]:dy, 0:z[-1]:dz] \n    # Reordering the axes and interpolating the image \n    new_grid = np.moveaxis(new_grid, (0, 1, 2, 3), (3, 0, 1, 2))\n    interpolated = f_scan(new_grid)\n    # We resize the interpolation to typical 512x512, and move the Z axis:\n    isotropic = np.zeros((interpolated.shape[2],512,512))\n    for k in range(interpolated.shape[2]):\n        interpolated_slice = interpolated[:,:,k]\n        interpolated_slice = cv2.resize(interpolated_slice, (512,512))\n        isotropic[k,:,:] = interpolated_slice\n\n    return(isotropic) \n\n\n# We open the CSV document and inspect the study names (which will be used to navigate throughout the RSPECT dataset and choose the indicated folders):\nDataFrame = pd.read_csv(\"../input/big-subset-pe/subset_RESPECT.csv\")\nfolders_of_interest = np.asarray(DataFrame)[:,1]\n\n# We also open the train CSV document (i.e., the original, this is intended in order to check labels and file names, as the previous CSV only includes \n# study-level information):\nFull_DataFrame = pd.read_csv(\"/kaggle/input/original-train-data/train_data_original_RSNA.csv\")\n\n# Creating folder for download (will include zip files):\nos.makedirs(\"./subset_zip\", exist_ok = True)\nos.makedirs(\"./GT_RVLV_zip\", exist_ok = True)\nos.makedirs(\"./GT_PE_zip\", exist_ok = True)\n\n# The train folder is full of several folders which, in turn, contain a folder with the patient's ID, and, inside the latter, all the different\n# independent .dcm files. We check each folder name in our list of interest and copy the folder it contains (i.e., the patient folder) in the \n# new directory:\ni = 0\nc = 0\nstart = time.time()\n\nfor folder in folders_of_interest:\n    if i >= (c*20) and i < ((c+1)*20) and i != 223 and i != 224 and i != 225 and i != 227:\n        GT_PE = []\n        path_to_folder = \"/kaggle/input/rsna-str-pulmonary-embolism-detection/train/\" + folder \n        # We enter the patient folder (named after the patient's ID) inside the general folder:\n        patient_folder_name = os.listdir(path_to_folder)\n        full_path_interest = path_to_folder + \"/\" + patient_folder_name[0]\n        # We get the list of files contained in the path:\n        files = os.listdir(full_path_interest)\n\n        # We create a new directory in the subset folder. The structure will be: a folder with the original folder name (i.e., first folder in train)\n        # with the files contained inside (we do not want this time a folder inside a folder):\n        if i < 10:\n            num_i = \"00\" + str(i)\n        elif 10 <= i and i < 100:\n            num_i = \"0\" + str(i)\n        else:\n            num_i = i\n\n        new_folder = \"./subset_zip/\" + str(num_i) + \"_\" + folder\n        os.makedirs(new_folder, exist_ok = True) \n\n        # We will build our GT lists (later on will get converted to numpy arrays):\n        count_slice = 0\n        for dcm_slice in files:\n            # Since we need to build our ground truth labels ordered, we will append them to a GT list that will be later on converted to an array. \n            # The labels of interest are: pe_present_on_image (4th feature in csv) and rv_lv_ratio_gte_1 (9th feature in csv).\n            info_slice = Full_DataFrame[Full_DataFrame.StudyInstanceUID == folder].iloc[count_slice]\n            count_slice = count_slice + 1\n\n            # Note: Rows in the CSV file are ordered while the actual images are disordered (Kaggle sorts them alphanumerically so both instances do \n            # not match). Later on, we will attend this problematic by renaming the images we read as we save them in a new (smaller) folder.\n            PE_present = info_slice[3]        \n            GT_PE.append(PE_present)\n\n        GT_PE = np.asarray(GT_PE)\n\n        # As I have noticed GT_PE includes values bigger than 1 (e.g., 2), maybe due to human errors, we will normalize that easily (and ensure 0 or 1 only).\n        GT_PE = 1*(GT_PE > 0)\n\n        namePE = str(num_i) + \"_GT_PE\"\n\n        RV_LV = info_slice[8] \n        GT_RVLV = np.asarray(RV_LV)\n        nameRVLV = str(num_i) + \"_GT_RVLV\"\n        np.save(\"/kaggle/working/GT_RVLV_zip/\" + nameRVLV, GT_RVLV)\n\n        # Now we will open, transform and save the diferent .dcm files as .png files:\n        # We will also create an array where the images will be introduced. The array \n        # will be transformed to guarantee isotropic voxels:\n        array_study = np.zeros((2000,512,512))\n        for dcm_slice in files:    \n            # We get each original file contained in train, convert it to a PNG image, and copy it in the new folder:\n            original_file_path = full_path_interest + \"/\" + dcm_slice\n            #img = sitk.ReadImage(original_file_path)\n            img = pydicom.read_file(original_file_path)\n            array = img.pixel_array\n\n            # We get the z-position of the slice from the metadata of the file (we will use this information to name the files and order the images):\n            ds = pydicom.read_file(original_file_path)\n            num_z = int(ds.InstanceNumber)\n\n            # We convert the raw data (array type) to Hounsfield Units (HU):\n            array = array*int(img.RescaleSlope) + np.ones((array.shape))*int(img.RescaleIntercept)\n            #array = sitk.GetArrayFromImage(img)\n            #array = array[0,:,:]\n\n            # Windowing (clipping):\n            lower_HU = -250\n            upper_HU = 450        \n            subtraction = np.subtract(array,(lower_HU*np.ones(array.shape)))\n            array = (subtraction/(upper_HU-lower_HU))    \n            array[array <= 0] = 0\n            array[array >= 1] = 1\n            array = array*255\n\n            # Resize:\n            array = cv2.resize(array, (512,512))\n\n            # Convert 16-bit pixels to 8-bit:\n            array = array.astype(np.uint8)\n\n            # Introducing the array into the general array:\n            array_study[num_z-1,:,:] = array\n\n        # We get only the elements of the array_study that contain information. Note that, in this big array, the elements are ordered:\n        aux_small_array = np.zeros((len(GT_PE),512,512))\n        for n in range(len(array_study)):\n            if np.sum(array_study[n,:,:]) != 0:\n                low_limit = n\n                break\n        aux_small_array = array_study[n:(n + len(GT_PE)),:,:]\n        array_study = aux_small_array\n\n        # Isotropic voxels:\n        # Original DICOM thickness and XY spacing:\n        thick = int(ds.SliceThickness)\n        spaceX = float(ds.PixelSpacing[0])\n        spaceY = float(ds.PixelSpacing[1])\n\n        # Transformation to isotropy:\n        array_study = do_interpolate(array_study, [spaceX, spaceY, thick])\n\n        # We also transform the PE array so it can match the new array:                            \n        old_x = np.linspace(0, len(GT_PE), num = len(GT_PE))    \n        new_x = np.linspace(0, len(GT_PE), num = len(array_study))\n\n        f = interpolate.interp1d(old_x, GT_PE, kind = \"nearest-up\")\n        resampled_GT_PE = f(new_x)\n        resampled_GT_PE = resampled_GT_PE.astype(np.uint8)\n\n        # Saving the PE_present labels (new):\n        np.save(\"/kaggle/working/GT_PE_zip/\" + namePE, resampled_GT_PE)\n\n        # Getting proper name for images and saving each PNG independently:\n        for array_index in range(len(array_study)):\n            array = array_study[array_index,:,:]\n            if array_index < 10:\n                name = \"00\" + str(array_index)\n            elif 10 <= array_index and array_index < 100:\n                name = \"0\" + str(array_index)\n            else:\n                name = str(array_index)\n            final_image_path = new_folder + \"/\" + name + \".png\"\n            cv2.imwrite(final_image_path, array)\n        print(\"Iteration number \", i, \" complete.\")\n        end = time.time()\n        print(end - start)\n    i = i + 1\n    \n    #if i == 10:  # Pause to get only 10 studies (micro_set for small checks).\n     #   break\n       \n        \n# We create a zip file with the subset folder we will be using later on:\nmake_archive(base_name = 'download_subset_png', format = 'zip', root_dir = \"./subset_zip\")\n\n# We also create zips for the arrays:\nmake_archive(base_name = 'download_zip_RVLV', format = 'zip', root_dir = \"./GT_RVLV_zip\")\nmake_archive(base_name = 'download_zip_PE', format = 'zip', root_dir = \"./GT_PE_zip\")\n\nimport IPython.display as ipd\nbeep = np.sin(2*np.pi*400*np.arange(10000*2)/10000)\nipd.Audio(beep, rate=10000, autoplay=True)\n'''","metadata":{"execution":{"iopub.status.busy":"2023-05-21T08:22:59.511555Z","iopub.execute_input":"2023-05-21T08:22:59.51215Z","iopub.status.idle":"2023-05-21T08:25:17.271868Z","shell.execute_reply.started":"2023-05-21T08:22:59.512067Z","shell.execute_reply":"2023-05-21T08:25:17.26991Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Relevant imports:\nimport scipy\nimport scipy.interpolate as interpolate\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nfrom shutil import copy, make_archive\nimport pydicom\nimport matplotlib.pyplot as plt\nimport cv2\nimport SimpleITK as sitk\nimport os\nimport copy\nimport time\n\n# We create a function that takes a CT and makes its voxels isotropic:\ndef do_interpolate(image_data, steps):\n    dx, dy, dz = 1., 1., 1.   # Desired (isotropic) spacing\n    image_data = np.moveaxis(image_data, 0, -1)\n    x, y, z = [steps[k] * np.arange(image_data.shape[k]) for k in range(len(steps))]  # original grid\n    # We create an interpolator:\n    f_scan = interpolate.RegularGridInterpolator((x, y, z), image_data, bounds_error = False, method = \"nearest\") \n    # New grid (with adequate spacing):   \n    new_grid = np.mgrid[0:x[-1]:dx, 0:y[-1]:dy, 0:z[-1]:dz] \n    # Reordering the axes and interpolating the image \n    new_grid = np.moveaxis(new_grid, (0, 1, 2, 3), (3, 0, 1, 2))\n    interpolated = f_scan(new_grid)\n    # We resize the interpolation to typical 512x512, and move the Z axis:\n    isotropic = np.zeros((interpolated.shape[2],512,512))\n    for k in range(interpolated.shape[2]):\n        interpolated_slice = interpolated[:,:,k]\n        interpolated_slice = cv2.resize(interpolated_slice, (512,512))\n        isotropic[k,:,:] = interpolated_slice\n\n    return(isotropic) \n\n# We open the CSV document and inspect the study names (which will be used to navigate throughout the RSPECT dataset and choose the indicated folders):\nDataFrame = pd.read_csv(\"../input/test-subset/subset_RESPECT.csv\")\nfolders_of_interest = np.asarray(DataFrame)[:,1]\n\n# We also open the train CSV document (i.e., the original, this is intended in order to check labels and file names, as the previous CSV only includes \n# study-level information):\nFull_DataFrame = pd.read_csv(\"/kaggle/input/original-train-data/train_data_original_RSNA.csv\")\n\n# Creating folder for download (will include zip files):\nos.makedirs(\"./subset_zip\", exist_ok = True)\nos.makedirs(\"./GT_RVLV_zip\", exist_ok = True)\nos.makedirs(\"./GT_PE_zip\", exist_ok = True)\n\n# The train folder is full of several folders which, in turn, contain a folder with the patient's ID, and, inside the latter, all the different\n# independent .dcm files. We check each folder name in our list of interest and copy the folder it contains (i.e., the patient folder) in the \n# new directory:\ni = 0\nc = 0\nfor folder in folders_of_interest:\n    if i >= c*100 and i < (c+1)*100 and i != 195 and i != 214 and i != 223 and i != 224 and i != 225 and i != 227 and i != 459 and i != 470 and i != 546 and i != 674 and i != 684:\n        #print(folder)\n        path_to_folder = \"/kaggle/input/rsna-str-pulmonary-embolism-detection/train/\" + folder \n        # We enter the patient folder (named after the patient's ID) inside the general folder:\n        patient_folder_name = os.listdir(path_to_folder)\n        full_path_interest = path_to_folder + \"/\" + patient_folder_name[0]\n        #print(patient_folder_name[0])\n        # We get the list of files contained in the path:\n        files = os.listdir(full_path_interest)\n\n        GT_PE = np.zeros(len(files))\n\n        # Getting minimum:\n        min_z = 100000\n        for dcm_slice in files:    \n            # We get each original file contained in train, convert it to a PNG image, and copy it in the new folder:\n            original_file_path = full_path_interest + \"/\" + dcm_slice\n            #img = sitk.ReadImage(original_file_path)\n            img = pydicom.read_file(original_file_path)\n            array = img.pixel_array\n\n            # We get the z-position of the slice from the metadata of the file (we will use this information to name the files and order the images):\n            ds = pydicom.read_file(original_file_path)\n            num_z = int(ds.InstanceNumber)\n            if num_z - 1 <= min_z:\n                min_z = copy.copy(num_z - 1)\n       \n\n        # We create a new directory in the subset folder. The structure will be: a folder with the original folder name (i.e., first folder in train)\n        # with the files contained inside (we do not want this time a folder inside a folder):\n        if i < 10:\n            num_i = \"00\" + str(i)\n        elif 10 <= i and i < 100:\n            num_i = \"0\" + str(i)\n        else:\n            num_i = i\n\n        new_folder = \"./subset_zip/\" + str(num_i) + \"_\" + folder\n        os.makedirs(new_folder, exist_ok = True) \n\n        # We will build our GT lists (later on will get converted to numpy arrays):\n        count_slice = 0\n        for dcm_slice in files:\n            dcm_slice = dcm_slice.replace(\".dcm\",\"\")\n            # Since we need to build our ground truth labels ordered, we will append them to a GT list that will be later on converted to an array. \n            # The labels of interest are: pe_present_on_image (4th feature in csv) and rv_lv_ratio_gte_1 (9th feature in csv).\n            info_slice = Full_DataFrame[Full_DataFrame.SOPInstanceUID == dcm_slice]\n            count_slice = count_slice + 1\n\n            # Note: Rows in the CSV file are ordered while the actual images are disordered (Kaggle sorts them alphanumerically so both instances do \n            # not match). Later on, we will attend this problematic by renaming the images we read as we save them in a new (smaller) folder.\n            PE_present = np.asarray(info_slice)[0][3]        \n\n            # Position of slice:\n            original_file_path = full_path_interest + \"/\" + dcm_slice + \".dcm\"\n            ds = pydicom.read_file(original_file_path)\n            num_z = int(ds.InstanceNumber)\n\n            # Including PE_label in the correct position:\n            GT_PE[(num_z-1)-min_z] = PE_present\n\n        #GT_PE = np.asarray(GT_PE)\n\n        RV_LV = np.asarray(info_slice)[0][8] \n        GT_RVLV = np.asarray(RV_LV)\n        nameRVLV = str(num_i) + \"_GT_RVLV\"\n        np.save(\"/kaggle/working/GT_RVLV_zip/\" + nameRVLV, GT_RVLV)\n        \n        # As I have noticed GT_PE includes values bigger than 1 (e.g., 2), maybe due to human errors, we will normalize that easily (and ensure 0 or 1 only).\n        GT_PE = 1*(GT_PE > 0)\n\n        namePE = str(num_i) + \"_GT_PE\"\n\n        # Now we will open, transform and save the diferent .dcm files as .png files:\n        # We will also create an array where the images will be introduced. The array \n        # will be transformed to guarantee isotropic voxels:\n        array_study = np.zeros((2000,512,512))\n        for dcm_slice in files:    \n            # We get each original file contained in train, convert it to a PNG image, and copy it in the new folder:\n            original_file_path = full_path_interest + \"/\" + dcm_slice\n            #img = sitk.ReadImage(original_file_path)\n            img = pydicom.read_file(original_file_path)\n            array = img.pixel_array\n\n            # We get the z-position of the slice from the metadata of the file (we will use this information to name the files and order the images):\n            ds = pydicom.read_file(original_file_path)\n            num_z = int(ds.InstanceNumber)\n\n            # We convert the raw data (array type) to Hounsfield Units (HU):\n            array = array*int(img.RescaleSlope) + np.ones((array.shape))*int(img.RescaleIntercept)\n            #array = sitk.GetArrayFromImage(img)\n            #array = array[0,:,:]\n\n            # Windowing (clipping):\n            lower_HU = -250\n            upper_HU = 450        \n            subtraction = np.subtract(array,(lower_HU*np.ones(array.shape)))\n            array = (subtraction/(upper_HU-lower_HU))    \n            array[array <= 0] = 0\n            array[array >= 1] = 1\n            array = array*255\n\n            # Resize:\n            array = cv2.resize(array, (512,512))\n\n            # Convert 16-bit pixels to 8-bit:\n            array = array.astype(np.uint8)\n\n            # Introducing the array into the general array:\n            array_study[num_z-1,:,:] = array\n\n        # We get only the elements of the array_study that contain information. Note that, in this big array, the elements are ordered:\n        aux_small_array = np.zeros((len(GT_PE),512,512))\n        for n in range(len(array_study)):\n            if np.sum(array_study[n,:,:]) != 0:\n                low_limit = n\n                break\n        aux_small_array = array_study[n:(n + len(GT_PE)),:,:]\n        array_study = aux_small_array\n\n        # Isotropic voxels:\n        # Original DICOM thickness and XY spacing:\n        thick = int(ds.SliceThickness)\n        spaceX = float(ds.PixelSpacing[0])\n        spaceY = float(ds.PixelSpacing[1])\n\n        # Transformation to isotropy:\n        array_study = do_interpolate(array_study, [spaceX, spaceY, thick])\n\n        # We also transform the PE array so it can match the new array:                            \n        old_x = np.linspace(0, len(GT_PE), num = len(GT_PE))    \n        new_x = np.linspace(0, len(GT_PE), num = len(array_study))\n\n        f = interpolate.interp1d(old_x, GT_PE, kind = \"nearest-up\")\n        resampled_GT_PE = f(new_x)\n        resampled_GT_PE = resampled_GT_PE.astype(np.uint8)\n\n        # Saving the PE_present labels (new):\n        np.save(\"/kaggle/working/GT_PE_zip/\" + namePE, resampled_GT_PE)\n\n        # Getting proper name for images and saving each PNG independently:\n        for array_index in range(len(array_study)):\n            array = array_study[array_index,:,:]\n            if array_index < 10:\n                name = \"00\" + str(array_index)\n            elif 10 <= array_index and array_index < 100:\n                name = \"0\" + str(array_index)\n            else:\n                name = str(array_index)\n            final_image_path = new_folder + \"/\" + name + \".png\"\n            cv2.imwrite(final_image_path, array)\n        \n        print(\"Iteration number \", i, \" complete.\")\n    #break\n    i = i + 1\n    #if i == 5:  # Pause to get only 10 studies (micro_set for small checks).\n     #   break\n       \n        \n# We create a zip file with the subset folder we will be using later on:\n\n# We also create zips for the arrays:\nmake_archive(base_name = 'download_subset_png', format = 'zip', root_dir = \"./subset_zip\")\n\n# We also create zips for the arrays:\nmake_archive(base_name = 'download_zip_RVLV', format = 'zip', root_dir = \"./GT_RVLV_zip\")\nmake_archive(base_name = 'download_zip_PE', format = 'zip', root_dir = \"./GT_PE_zip\")","metadata":{"execution":{"iopub.status.busy":"2023-05-23T17:47:05.381712Z","iopub.execute_input":"2023-05-23T17:47:05.382238Z","iopub.status.idle":"2023-05-23T18:28:25.479796Z","shell.execute_reply.started":"2023-05-23T17:47:05.382205Z","shell.execute_reply":"2023-05-23T18:28:25.478613Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\n# 29/09/2023\n#-----------------------------------------------------------------------------------------\n\n# ~ GENERAL ~\n# Relevant imports:\nimport scipy\nimport scipy.interpolate as interpolate\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nfrom shutil import copy, make_archive\nimport pydicom\nimport matplotlib.pyplot as plt\nimport cv2\nimport SimpleITK as sitk\nimport os\nimport copy\nimport time\n\n# We create a function that takes a CT and makes its voxels isotropic:\ndef do_interpolate(image_data, steps):\n    dx, dy, dz = 1., 1., 1.   # Desired (isotropic) spacing\n    image_data = np.moveaxis(image_data, 0, -1)\n    x, y, z = [steps[k] * np.arange(image_data.shape[k]) for k in range(len(steps))]  # original grid\n    # We create an interpolator:\n    f_scan = interpolate.RegularGridInterpolator((x, y, z), image_data, bounds_error = False, method = \"nearest\") \n    # New grid (with adequate spacing):   \n    new_grid = np.mgrid[0:x[-1]:dx, 0:y[-1]:dy, 0:z[-1]:dz] \n    # Reordering the axes and interpolating the image \n    new_grid = np.moveaxis(new_grid, (0, 1, 2, 3), (3, 0, 1, 2))\n    interpolated = f_scan(new_grid)\n    # We resize the interpolation to typical 512x512, and move the Z axis:\n    isotropic = np.zeros((interpolated.shape[2],512,512))\n    for k in range(interpolated.shape[2]):\n        interpolated_slice = interpolated[:,:,k]\n        interpolated_slice = cv2.resize(interpolated_slice, (512,512))\n        isotropic[k,:,:] = interpolated_slice\n\n    return(isotropic) \n\n\n\n\n\n\n\n\n\n\n\n\n\n\n#-----------------------------------------------------------------------------------------\n# ~ TEST ~ \n# We open the CSV document and inspect the study names (which will be used to navigate throughout the RSPECT dataset and choose the indicated folders):\nDataFrame = pd.read_csv(\"../input/test-subset/test_subset_RESPECT.csv\")\nfolders_of_interest = np.asarray(DataFrame)[:,1]\n\n# We also open the train CSV document (i.e., the original, this is intended in order to check labels and file names, as the previous CSV only includes \n# study-level information):\nFull_DataFrame = pd.read_csv(\"/kaggle/input/original-train-data/train_data_original_RSNA.csv\")\n\n# Creating folder for download (will include zip files):\nos.makedirs(\"./subset_zip\", exist_ok = True)\nos.makedirs(\"./GT_RVLV_zip\", exist_ok = True)\nos.makedirs(\"./GT_PE_zip\", exist_ok = True)\n\n# The train folder is full of several folders which, in turn, contain a folder with the patient's ID, and, inside the latter, all the different\n# independent .dcm files. We check each folder name in our list of interest and copy the folder it contains (i.e., the patient folder) in the \n# new directory:\ni = 0\nc = 0\nfor folder in folders_of_interest:\n    if i >= c*100 and i < (c+1)*100 and i != 195 and i != 214 and i != 223 and i != 224 and i != 225 and i != 227 and i != 459 and i != 470 and i != 546 and i != 674 and i != 684:\n        #print(folder)\n        path_to_folder = \"/kaggle/input/rsna-str-pulmonary-embolism-detection/train/\" + folder \n        # We enter the patient folder (named after the patient's ID) inside the general folder:\n        patient_folder_name = os.listdir(path_to_folder)\n        full_path_interest = path_to_folder + \"/\" + patient_folder_name[0]\n        #print(patient_folder_name[0])\n        # We get the list of files contained in the path:\n        files = os.listdir(full_path_interest)\n\n        GT_PE = np.zeros(len(files))\n\n        # Getting minimum:\n        min_z = 100000\n        for dcm_slice in files:    \n            # We get each original file contained in train, convert it to a PNG image, and copy it in the new folder:\n            original_file_path = full_path_interest + \"/\" + dcm_slice\n            #img = sitk.ReadImage(original_file_path)\n            img = pydicom.read_file(original_file_path)\n            array = img.pixel_array\n\n            # We get the z-position of the slice from the metadata of the file (we will use this information to name the files and order the images):\n            ds = pydicom.read_file(original_file_path)\n            num_z = int(ds.InstanceNumber)\n            if num_z - 1 <= min_z:\n                min_z = copy.copy(num_z - 1)\n       \n\n        # We create a new directory in the subset folder. The structure will be: a folder with the original folder name (i.e., first folder in train)\n        # with the files contained inside (we do not want this time a folder inside a folder):\n        if i < 10:\n            num_i = \"00\" + str(i)\n        elif 10 <= i and i < 100:\n            num_i = \"0\" + str(i)\n        else:\n            num_i = i\n\n        new_folder = \"./subset_zip/\" + str(num_i) + \"_\" + folder\n        os.makedirs(new_folder, exist_ok = True) \n\n        # We will build our GT lists (later on will get converted to numpy arrays):\n        count_slice = 0\n        for dcm_slice in files:\n            dcm_slice = dcm_slice.replace(\".dcm\",\"\")\n            # Since we need to build our ground truth labels ordered, we will append them to a GT list that will be later on converted to an array. \n            # The labels of interest are: pe_present_on_image (4th feature in csv) and rv_lv_ratio_gte_1 (9th feature in csv).\n            info_slice = Full_DataFrame[Full_DataFrame.SOPInstanceUID == dcm_slice]\n            count_slice = count_slice + 1\n\n            # Note: Rows in the CSV file are ordered while the actual images are disordered (Kaggle sorts them alphanumerically so both instances do \n            # not match). Later on, we will attend this problematic by renaming the images we read as we save them in a new (smaller) folder.\n            PE_present = np.asarray(info_slice)[0][3]        \n\n            # Position of slice:\n            original_file_path = full_path_interest + \"/\" + dcm_slice + \".dcm\"\n            ds = pydicom.read_file(original_file_path)\n            num_z = int(ds.InstanceNumber)\n\n            # Including PE_label in the correct position:\n            GT_PE[(num_z-1)-min_z] = PE_present\n\n        #GT_PE = np.asarray(GT_PE)\n\n        RV_LV = np.asarray(info_slice)[0][8] \n        GT_RVLV = np.asarray(RV_LV)\n        nameRVLV = str(num_i) + \"_GT_RVLV\"\n        np.save(\"/kaggle/working/GT_RVLV_zip/\" + nameRVLV, GT_RVLV)\n        \n        # As I have noticed GT_PE includes values bigger than 1 (e.g., 2), maybe due to human errors, we will normalize that easily (and ensure 0 or 1 only).\n        GT_PE = 1*(GT_PE > 0)\n\n        namePE = str(num_i) + \"_GT_PE\"\n\n        # Now we will open, transform and save the diferent .dcm files as .png files:\n        # We will also create an array where the images will be introduced. The array \n        # will be transformed to guarantee isotropic voxels:\n        array_study = np.zeros((2000,512,512))\n        for dcm_slice in files:    \n            # We get each original file contained in train, convert it to a PNG image, and copy it in the new folder:\n            original_file_path = full_path_interest + \"/\" + dcm_slice\n            #img = sitk.ReadImage(original_file_path)\n            img = pydicom.read_file(original_file_path)\n            array = img.pixel_array\n\n            # We get the z-position of the slice from the metadata of the file (we will use this information to name the files and order the images):\n            ds = pydicom.read_file(original_file_path)\n            num_z = int(ds.InstanceNumber)\n\n            # We convert the raw data (array type) to Hounsfield Units (HU):\n            array = array*int(img.RescaleSlope) + np.ones((array.shape))*int(img.RescaleIntercept)\n            #array = sitk.GetArrayFromImage(img)\n            #array = array[0,:,:]\n\n            # Windowing (clipping):\n            lower_HU = -250\n            upper_HU = 450        \n            subtraction = np.subtract(array,(lower_HU*np.ones(array.shape)))\n            array = (subtraction/(upper_HU-lower_HU))    \n            array[array <= 0] = 0\n            array[array >= 1] = 1\n            array = array*255\n\n            # Resize:\n            array = cv2.resize(array, (512,512))\n\n            # Convert 16-bit pixels to 8-bit:\n            array = array.astype(np.uint8)\n\n            # Introducing the array into the general array:\n            array_study[num_z-1,:,:] = array\n\n        # We get only the elements of the array_study that contain information. Note that, in this big array, the elements are ordered:\n        aux_small_array = np.zeros((len(GT_PE),512,512))\n        for n in range(len(array_study)):\n            if np.sum(array_study[n,:,:]) != 0:\n                low_limit = n\n                break\n        aux_small_array = array_study[n:(n + len(GT_PE)),:,:]\n        array_study = aux_small_array\n\n        # Isotropic voxels:\n        # Original DICOM thickness and XY spacing:\n        thick = int(ds.SliceThickness)\n        spaceX = float(ds.PixelSpacing[0])\n        spaceY = float(ds.PixelSpacing[1])\n\n        # Transformation to isotropy:\n        array_study = do_interpolate(array_study, [spaceX, spaceY, thick])\n\n        # We also transform the PE array so it can match the new array:                            \n        old_x = np.linspace(0, len(GT_PE), num = len(GT_PE))    \n        new_x = np.linspace(0, len(GT_PE), num = len(array_study))\n\n        f = interpolate.interp1d(old_x, GT_PE, kind = \"nearest-up\")\n        resampled_GT_PE = f(new_x)\n        resampled_GT_PE = resampled_GT_PE.astype(np.uint8)\n\n        # Saving the PE_present labels (new):\n        np.save(\"/kaggle/working/GT_PE_zip/\" + namePE, resampled_GT_PE)\n\n        # Getting proper name for images and saving each PNG independently:\n        for array_index in range(len(array_study)):\n            array = array_study[array_index,:,:]\n            if array_index < 10:\n                name = \"00\" + str(array_index)\n            elif 10 <= array_index and array_index < 100:\n                name = \"0\" + str(array_index)\n            else:\n                name = str(array_index)\n            final_image_path = new_folder + \"/\" + name + \".png\"\n            cv2.imwrite(final_image_path, array)\n        \n        print(\"Iteration number \", i, \" complete.\")\n    #break\n    i = i + 1\n    #if i == 5:  # Pause to get only 10 studies (micro_set for small checks).\n     #   break\n       \n        \n# We create a zip file with the subset folder we will be using later on:\n\n# We also create zips for the arrays:\nmake_archive(base_name = 'download_subset_png', format = 'zip', root_dir = \"./subset_zip\")\n\n# We also create zips for the arrays:\nmake_archive(base_name = 'download_zip_RVLV', format = 'zip', root_dir = \"./GT_RVLV_zip\")\nmake_archive(base_name = 'download_zip_PE', format = 'zip', root_dir = \"./GT_PE_zip\")\n\n\n\n\n\n\n\n\n\n#------------------------------------------------------------------------------------------\n# ~ TRAIN ~ \n# We open the CSV document and inspect the study names (which will be used to navigate throughout the RSPECT dataset and choose the indicated folders):\nDataFrame = pd.read_csv(\"/kaggle/input/big-subset-pe\")\nfolders_of_interest = np.asarray(DataFrame)[:,1]\n\n# We also open the train CSV document (i.e., the original, this is intended in order to check labels and file names, as the previous CSV only includes \n# study-level information):\nFull_DataFrame = pd.read_csv(\"/kaggle/input/original-train-data/train_data_original_RSNA.csv\")\n\n# Creating folder for download (will include zip files):\nos.makedirs(\"./subset_zip_TRAIN\", exist_ok = True)\nos.makedirs(\"./GT_RVLV_zip_TRAIN\", exist_ok = True)\nos.makedirs(\"./GT_PE_zip_TRAIN\", exist_ok = True)\n\n# The train folder is full of several folders which, in turn, contain a folder with the patient's ID, and, inside the latter, all the different\n# independent .dcm files. We check each folder name in our list of interest and copy the folder it contains (i.e., the patient folder) in the \n# new directory:\ni = 0\nc = 0\nfor folder in folders_of_interest:\n    if i >= c*100 and i < (c+1)*100 and i != 195 and i != 214 and i != 223 and i != 224 and i != 225 and i != 227 and i != 459 and i != 470 and i != 546 and i != 674 and i != 684:\n        #print(folder)\n        path_to_folder = \"/kaggle/input/rsna-str-pulmonary-embolism-detection/train/\" + folder \n        # We enter the patient folder (named after the patient's ID) inside the general folder:\n        patient_folder_name = os.listdir(path_to_folder)\n        full_path_interest = path_to_folder + \"/\" + patient_folder_name[0]\n        #print(patient_folder_name[0])\n        # We get the list of files contained in the path:\n        files = os.listdir(full_path_interest)\n\n        GT_PE = np.zeros(len(files))\n\n        # Getting minimum:\n        min_z = 100000\n        for dcm_slice in files:    \n            # We get each original file contained in train, convert it to a PNG image, and copy it in the new folder:\n            original_file_path = full_path_interest + \"/\" + dcm_slice\n            #img = sitk.ReadImage(original_file_path)\n            img = pydicom.read_file(original_file_path)\n            array = img.pixel_array\n\n            # We get the z-position of the slice from the metadata of the file (we will use this information to name the files and order the images):\n            ds = pydicom.read_file(original_file_path)\n            num_z = int(ds.InstanceNumber)\n            if num_z - 1 <= min_z:\n                min_z = copy.copy(num_z - 1)\n       \n\n        # We create a new directory in the subset folder. The structure will be: a folder with the original folder name (i.e., first folder in train)\n        # with the files contained inside (we do not want this time a folder inside a folder):\n        if i < 10:\n            num_i = \"00\" + str(i)\n        elif 10 <= i and i < 100:\n            num_i = \"0\" + str(i)\n        else:\n            num_i = i\n\n        new_folder = \"./subset_zip/\" + str(num_i) + \"_\" + folder\n        os.makedirs(new_folder, exist_ok = True) \n\n        # We will build our GT lists (later on will get converted to numpy arrays):\n        count_slice = 0\n        for dcm_slice in files:\n            dcm_slice = dcm_slice.replace(\".dcm\",\"\")\n            # Since we need to build our ground truth labels ordered, we will append them to a GT list that will be later on converted to an array. \n            # The labels of interest are: pe_present_on_image (4th feature in csv) and rv_lv_ratio_gte_1 (9th feature in csv).\n            info_slice = Full_DataFrame[Full_DataFrame.SOPInstanceUID == dcm_slice]\n            count_slice = count_slice + 1\n\n            # Note: Rows in the CSV file are ordered while the actual images are disordered (Kaggle sorts them alphanumerically so both instances do \n            # not match). Later on, we will attend this problematic by renaming the images we read as we save them in a new (smaller) folder.\n            PE_present = np.asarray(info_slice)[0][3]        \n\n            # Position of slice:\n            original_file_path = full_path_interest + \"/\" + dcm_slice + \".dcm\"\n            ds = pydicom.read_file(original_file_path)\n            num_z = int(ds.InstanceNumber)\n\n            # Including PE_label in the correct position:\n            GT_PE[(num_z-1)-min_z] = PE_present\n\n        #GT_PE = np.asarray(GT_PE)\n\n        RV_LV = np.asarray(info_slice)[0][8] \n        GT_RVLV = np.asarray(RV_LV)\n        nameRVLV = str(num_i) + \"_GT_RVLV\"\n        np.save(\"/kaggle/working/GT_RVLV_zip/\" + nameRVLV, GT_RVLV)\n        \n        # As I have noticed GT_PE includes values bigger than 1 (e.g., 2), maybe due to human errors, we will normalize that easily (and ensure 0 or 1 only).\n        GT_PE = 1*(GT_PE > 0)\n\n        namePE = str(num_i) + \"_GT_PE\"\n\n        # Now we will open, transform and save the diferent .dcm files as .png files:\n        # We will also create an array where the images will be introduced. The array \n        # will be transformed to guarantee isotropic voxels:\n        array_study = np.zeros((2000,512,512))\n        for dcm_slice in files:    \n            # We get each original file contained in train, convert it to a PNG image, and copy it in the new folder:\n            original_file_path = full_path_interest + \"/\" + dcm_slice\n            #img = sitk.ReadImage(original_file_path)\n            img = pydicom.read_file(original_file_path)\n            array = img.pixel_array\n\n            # We get the z-position of the slice from the metadata of the file (we will use this information to name the files and order the images):\n            ds = pydicom.read_file(original_file_path)\n            num_z = int(ds.InstanceNumber)\n\n            # We convert the raw data (array type) to Hounsfield Units (HU):\n            array = array*int(img.RescaleSlope) + np.ones((array.shape))*int(img.RescaleIntercept)\n            #array = sitk.GetArrayFromImage(img)\n            #array = array[0,:,:]\n\n            # Windowing (clipping):\n            lower_HU = -250\n            upper_HU = 450        \n            subtraction = np.subtract(array,(lower_HU*np.ones(array.shape)))\n            array = (subtraction/(upper_HU-lower_HU))    \n            array[array <= 0] = 0\n            array[array >= 1] = 1\n            array = array*255\n\n            # Resize:\n            array = cv2.resize(array, (512,512))\n\n            # Convert 16-bit pixels to 8-bit:\n            array = array.astype(np.uint8)\n\n            # Introducing the array into the general array:\n            array_study[num_z-1,:,:] = array\n\n        # We get only the elements of the array_study that contain information. Note that, in this big array, the elements are ordered:\n        aux_small_array = np.zeros((len(GT_PE),512,512))\n        for n in range(len(array_study)):\n            if np.sum(array_study[n,:,:]) != 0:\n                low_limit = n\n                break\n        aux_small_array = array_study[n:(n + len(GT_PE)),:,:]\n        array_study = aux_small_array\n\n        # Isotropic voxels:\n        # Original DICOM thickness and XY spacing:\n        thick = int(ds.SliceThickness)\n        spaceX = float(ds.PixelSpacing[0])\n        spaceY = float(ds.PixelSpacing[1])\n\n        # Transformation to isotropy:\n        array_study = do_interpolate(array_study, [spaceX, spaceY, thick])\n\n        # We also transform the PE array so it can match the new array:                            \n        old_x = np.linspace(0, len(GT_PE), num = len(GT_PE))    \n        new_x = np.linspace(0, len(GT_PE), num = len(array_study))\n\n        f = interpolate.interp1d(old_x, GT_PE, kind = \"nearest-up\")\n        resampled_GT_PE = f(new_x)\n        resampled_GT_PE = resampled_GT_PE.astype(np.uint8)\n\n        # Saving the PE_present labels (new):\n        np.save(\"/kaggle/working/GT_PE_zip/\" + namePE, resampled_GT_PE)\n\n        # Getting proper name for images and saving each PNG independently:\n        for array_index in range(len(array_study)):\n            array = array_study[array_index,:,:]\n            if array_index < 10:\n                name = \"00\" + str(array_index)\n            elif 10 <= array_index and array_index < 100:\n                name = \"0\" + str(array_index)\n            else:\n                name = str(array_index)\n            final_image_path = new_folder + \"/\" + name + \".png\"\n            cv2.imwrite(final_image_path, array)\n        \n        print(\"Iteration number \", i, \" complete.\")\n    #break\n    i = i + 1\n    #if i == 5:  # Pause to get only 10 studies (micro_set for small checks).\n     #   break\n       \n        \n# We create a zip file with the subset folder we will be using later on:\n\n# We also create zips for the arrays:\nmake_archive(base_name = 'download_subset_png', format = 'zip', root_dir = \"./subset_zip\")\n\n# We also create zips for the arrays:\nmake_archive(base_name = 'download_zip_RVLV', format = 'zip', root_dir = \"./GT_RVLV_zip\")\nmake_archive(base_name = 'download_zip_PE', format = 'zip', root_dir = \"./GT_PE_zip\")\n'''","metadata":{},"execution_count":null,"outputs":[]}]}