{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"name":"Untitled0.ipynb","provenance":[],"collapsed_sections":[],"authorship_tag":"ABX9TyNof3khRW7RpbtyRO4LLdKk"},"kernelspec":{"name":"python3","display_name":"Python 3"}},"cells":[{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"Y5WzksAhSD26","executionInfo":{"status":"ok","timestamp":1612251630781,"user_tz":300,"elapsed":12952,"user":{"displayName":"Matthew So","photoUrl":"","userId":"18273672785190648211"}},"outputId":"79413b17-23c0-4cc9-857d-b7fbdea57e04"},"source":["#!pip install pytorch_lightning\r\n","#!pip install torchsummaryX\r\n","!pip install webdataset\r\n","# !pip install datasets\r\n","# !pip install wandb\r\n","#!pip install -r MedicalZooPytorch/installation/requirements.txt\r\n","#!pip install torch==1.7.1+cu101 torchvision==0.8.2+cu101 torchaudio==0.7.2 -f https://download.pytorch.org/whl/torch_stable.html\r\n","!git clone https://github.com/McMasterAI/Radiology-and-AI.git   #--branch augmentation \r\n","!git clone https://github.com/jcreinhold/intensity-normalization.git\r\n","! python intensity-normalization/setup.py install\r\n","!pip install scikit-fuzzy"],"execution_count":null,"outputs":[{"output_type":"stream","text":["Collecting webdataset\n","  Downloading https://files.pythonhosted.org/packages/16/bc/8b98d07eb97a51584ff305b468a13628f3905964d597c62513f6beacb4a4/webdataset-0.1.40-py3-none-any.whl\n","Collecting braceexpand\n","  Downloading https://files.pythonhosted.org/packages/bc/63/e944b82d7dc4fbba0199af90f6712626224990fbe7fd3fe323d9a984a679/braceexpand-0.1.6-py2.py3-none-any.whl\n","Requirement already satisfied: msgpack in /usr/local/lib/python3.6/dist-packages (from webdataset) (1.0.2)\n","Requirement already satisfied: torch in /usr/local/lib/python3.6/dist-packages (from webdataset) (1.7.0+cu101)\n","Collecting simplejson\n","\u001b[?25l  Downloading https://files.pythonhosted.org/packages/73/96/1e6b19045375890068d7342cbe280dd64ae73fd90b9735b5efb8d1e044a1/simplejson-3.17.2-cp36-cp36m-manylinux2010_x86_64.whl (127kB)\n","\u001b[K     |████████████████████████████████| 133kB 14.6MB/s \n","\u001b[?25hCollecting objectio\n","  Downloading https://files.pythonhosted.org/packages/86/e3/a132a91c4e9fd5e59c947263c7ef4e3415640fa151344f858e2def8c1726/objectio-0.2.29-py3-none-any.whl\n","Requirement already satisfied: Pillow in /usr/local/lib/python3.6/dist-packages (from webdataset) (7.0.0)\n","Requirement already satisfied: numpy in /usr/local/lib/python3.6/dist-packages (from webdataset) (1.19.5)\n","Requirement already satisfied: pyyaml in /usr/local/lib/python3.6/dist-packages (from webdataset) (3.13)\n","Requirement already satisfied: future in /usr/local/lib/python3.6/dist-packages (from torch->webdataset) (0.16.0)\n","Requirement already satisfied: typing-extensions in /usr/local/lib/python3.6/dist-packages (from torch->webdataset) (3.7.4.3)\n","Requirement already satisfied: dataclasses in /usr/local/lib/python3.6/dist-packages (from torch->webdataset) (0.8)\n","Collecting typer\n","  Downloading https://files.pythonhosted.org/packages/90/34/d138832f6945432c638f32137e6c79a3b682f06a63c488dcfaca6b166c64/typer-0.3.2-py3-none-any.whl\n","Requirement already satisfied: click<7.2.0,>=7.1.1 in /usr/local/lib/python3.6/dist-packages (from typer->objectio->webdataset) (7.1.2)\n","Installing collected packages: braceexpand, simplejson, typer, objectio, webdataset\n","Successfully installed braceexpand-0.1.6 objectio-0.2.29 simplejson-3.17.2 typer-0.3.2 webdataset-0.1.40\n","Cloning into 'Radiology-and-AI'...\n","remote: Enumerating objects: 152, done.\u001b[K\n","remote: Counting objects: 100% (152/152), done.\u001b[K\n","remote: Compressing objects: 100% (98/98), done.\u001b[K\n","remote: Total 312 (delta 92), reused 98 (delta 44), pack-reused 160\u001b[K\n","Receiving objects: 100% (312/312), 25.81 MiB | 23.06 MiB/s, done.\n","Resolving deltas: 100% (170/170), done.\n","Cloning into 'intensity-normalization'...\n","remote: Enumerating objects: 184, done.\u001b[K\n","remote: Counting objects: 100% (184/184), done.\u001b[K\n","remote: Compressing objects: 100% (167/167), done.\u001b[K\n","remote: Total 1216 (delta 111), reused 85 (delta 12), pack-reused 1032\u001b[K\n","Receiving objects: 100% (1216/1216), 1.03 MiB | 3.33 MiB/s, done.\n","Resolving deltas: 100% (846/846), done.\n","Traceback (most recent call last):\n","  File \"intensity-normalization/setup.py\", line 39, in <module>\n","    with open('README.md', encoding=\"utf-8\") as f:\n","FileNotFoundError: [Errno 2] No such file or directory: 'README.md'\n","Collecting scikit-fuzzy\n","\u001b[?25l  Downloading https://files.pythonhosted.org/packages/6c/f0/5eb5dbe0fd8dfe7d4651a8f4e591a196623a22b9e5339101e559695b4f6c/scikit-fuzzy-0.4.2.tar.gz (993kB)\n","\u001b[K     |████████████████████████████████| 1.0MB 7.9MB/s \n","\u001b[?25hRequirement already satisfied: numpy>=1.6.0 in /usr/local/lib/python3.6/dist-packages (from scikit-fuzzy) (1.19.5)\n","Requirement already satisfied: scipy>=0.9.0 in /usr/local/lib/python3.6/dist-packages (from scikit-fuzzy) (1.4.1)\n","Requirement already satisfied: networkx>=1.9.0 in /usr/local/lib/python3.6/dist-packages (from scikit-fuzzy) (2.5)\n","Requirement already satisfied: decorator>=4.3.0 in /usr/local/lib/python3.6/dist-packages (from networkx>=1.9.0->scikit-fuzzy) (4.4.2)\n","Building wheels for collected packages: scikit-fuzzy\n","  Building wheel for scikit-fuzzy (setup.py) ... \u001b[?25l\u001b[?25hdone\n","  Created wheel for scikit-fuzzy: filename=scikit_fuzzy-0.4.2-cp36-none-any.whl size=894070 sha256=dc48fb5886148b425bce7c036cabbb78c2f23ba5e7a71fb4a0eb36607a18db2b\n","  Stored in directory: /root/.cache/pip/wheels/b9/4e/77/da79b16f64ef1738d95486e2731eea09d73e90a72465096600\n","Successfully built scikit-fuzzy\n","Installing collected packages: scikit-fuzzy\n","Successfully installed scikit-fuzzy-0.4.2\n"],"name":"stdout"}]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"yxb27aj9SGCK","executionInfo":{"status":"ok","timestamp":1612251710085,"user_tz":300,"elapsed":92240,"user":{"displayName":"Matthew So","photoUrl":"","userId":"18273672785190648211"}},"outputId":"60bc3a10-deb6-46d9-9b3e-f799889937e4"},"source":["from google.colab import drive\r\n","drive.mount('/content/drive', force_remount=True)"],"execution_count":null,"outputs":[{"output_type":"stream","text":["Mounted at /content/drive\n"],"name":"stdout"}]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"jvZdWCBgSG83","executionInfo":{"status":"ok","timestamp":1612251710086,"user_tz":300,"elapsed":92233,"user":{"displayName":"Matthew So","photoUrl":"","userId":"18273672785190648211"}},"outputId":"95869817-a207-4118-dbcb-0edfcbc608c9"},"source":["cd drive/MyDrive/MacAI"],"execution_count":null,"outputs":[{"output_type":"stream","text":["/content/drive/MyDrive/MacAI\n"],"name":"stdout"}]},{"cell_type":"code","metadata":{"id":"Q_9WRvmPSIAp"},"source":["import sys\r\n","sys.path.append('./Radiology-and-AI/Radiology_and_AI')\r\n","sys.path.append('./intensity-normalization')\r\n","import os\r\n","import torch\r\n","import numpy as np\r\n","import webdataset as wds"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"DCsLbOQWWJcg","executionInfo":{"status":"ok","timestamp":1612251721823,"user_tz":300,"elapsed":103956,"user":{"displayName":"Matthew So","photoUrl":"","userId":"18273672785190648211"}},"outputId":"26f7ee31-f3ac-4f90-d1a1-6d931f9d37d4"},"source":["import intensity_normalization"],"execution_count":null,"outputs":[{"output_type":"stream","text":["/usr/local/lib/python3.6/dist-packages/statsmodels/tools/_testing.py:19: FutureWarning: pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.\n","  import pandas.util.testing as tm\n"],"name":"stderr"}]},{"cell_type":"code","metadata":{"id":"4w__ktMcuIOH"},"source":["from io import BytesIO\r\n","from nibabel import FileHolder, Nifti1Image\r\n","import torch\r\n","import numpy as np\r\n","from scipy.interpolate import RegularGridInterpolator\r\n","from scipy.ndimage.filters import gaussian_filter\r\n","from time import time"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"5PD3dOm3t5tF"},"source":["import matplotlib.pyplot as plt\r\n","import seaborn as sns"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"wBz2v1Wlwxap"},"source":["from scipy.interpolate import interp1d"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"-5QnFmqFS258"},"source":["train_dataset = wds.Dataset(\"macai_datasets/brats/train/brats_train.tar.gz\")\r\n","eval_dataset = wds.Dataset(\"macai_datasets/brats/validation/brats_validation.tar.gz\")"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"17KznqwBraQ-"},"source":["def np_img_collator(batch):\r\n","    bytes_data_list = [list(batch[i].items())[1][1] for i in range(5)]\r\n","    bytes_data_keys = [list(batch[i].items())[0][1].split('_')[-1] for i in range(5)]\r\n","    bytes_data_dict = dict(zip(bytes_data_keys,bytes_data_list))\r\n","\r\n","    bb = BytesIO(bytes_data_dict['flair'])\r\n","    fh = FileHolder(fileobj=bb)\r\n","    f_flair = Nifti1Image.from_file_map({'header': fh, 'image':fh}).get_fdata()\r\n","    bb = BytesIO(bytes_data_dict['seg'])\r\n","    fh = FileHolder(fileobj=bb)\r\n","    f_seg = Nifti1Image.from_file_map({'header': fh, 'image':fh}).get_fdata()\r\n","    bb = BytesIO(bytes_data_dict['t1'])\r\n","    fh = FileHolder(fileobj=bb)\r\n","    f_t1 = Nifti1Image.from_file_map({'header': fh, 'image':fh}).get_fdata()\r\n","    bb = BytesIO(bytes_data_dict['t1ce'])\r\n","    fh = FileHolder(fileobj=bb)\r\n","    f_t1ce=Nifti1Image.from_file_map({'header':fh, 'image':fh}).get_fdata()\r\n","    bb = BytesIO(bytes_data_dict['t2'])\r\n","    fh = FileHolder(fileobj=bb)\r\n","    f_t2 =Nifti1Image.from_file_map({'header':fh, 'image':fh}).get_fdata()\r\n","\r\n","    padding = [(0, 0), (0, 0), (0, 0)]# last (2,3)\r\n","    f_flair = np.expand_dims(np.pad(f_flair, padding), axis=0)\r\n","    f_t1 = np.expand_dims(np.pad(f_t1, padding), axis=0)\r\n","    f_t2 = np.expand_dims(np.pad(f_t2, padding), axis=0)\r\n","    f_t1ce = np.expand_dims(np.pad(f_t1ce, padding), axis=0)\r\n","\r\n","    f_seg = np.pad(f_seg, padding)\r\n","    concat = np.concatenate([f_t1, f_t1ce, f_t2, f_flair], axis=0)\r\n","    f_seg = np.expand_dims(f_seg, axis=0)\r\n","\r\n","    return ([concat, f_seg])"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"2tc_4akNsCSv"},"source":["train_dataloader = torch.utils.data.DataLoader(train_dataset, batch_size=5,collate_fn=np_img_collator)"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"cmcngwI6rQAN"},"source":["def nyul_train_dataloader(dataloader, n_imgs = 4, i_min=1, i_max=99, i_s_min=1, i_s_max=100, l_percentile=10, u_percentile=90, step=20):\r\n","    \"\"\"\r\n","    determine the standard scale for the set of images\r\n","    Args:\r\n","        img_fns (list): set of NifTI MR image paths which are to be normalized\r\n","        mask_fns (list): set of corresponding masks (if not provided, estimated)\r\n","        i_min (float): minimum percentile to consider in the images\r\n","        i_max (float): maximum percentile to consider in the images\r\n","        i_s_min (float): minimum percentile on the standard scale\r\n","        i_s_max (float): maximum percentile on the standard scale\r\n","        l_percentile (int): middle percentile lower bound (e.g., for deciles 10)\r\n","        u_percentile (int): middle percentile upper bound (e.g., for deciles 90)\r\n","        step (int): step for middle percentiles (e.g., for deciles 10)\r\n","    Returns:\r\n","        standard_scale (np.ndarray): average landmark intensity for images\r\n","        percs (np.ndarray): array of all percentiles used\r\n","    \"\"\"\r\n","    percss = [np.concatenate(([i_min], np.arange(l_percentile, u_percentile+1, step), [i_max])) for _ in range(n_imgs)]\r\n","    standard_scales = [np.zeros(len(percss[0])) for _ in range(n_imgs)]\r\n","\r\n","    iteration = 1\r\n","    for all_img, seg_data in dataloader:\r\n","      print(iteration)\r\n","    \r\n","     # print(seg_data.shape)\r\n","      mask_data = seg_data\r\n","      mask_data[seg_data ==0] = 1\r\n","      mask_data = np.squeeze(mask_data, axis=0)\r\n","\r\n","      #mask_data[mask_data==2] = 0 # ignore edema\r\n","\r\n","      for i in range(n_imgs):\r\n","        img_data = all_img[i]\r\n","\r\n","        masked = img_data[mask_data > 0]\r\n","        landmarks = intensity_normalization.normalize.nyul.get_landmarks(masked, percss[i])\r\n","        min_p = np.percentile(masked, i_min)\r\n","        max_p = np.percentile(masked, i_max)\r\n","        f = interp1d([min_p, max_p], [i_s_min, i_s_max])\r\n","        landmarks = np.array(f(landmarks))\r\n","        \r\n","        standard_scales[i] += landmarks\r\n","      iteration += 1\r\n","\r\n","    standard_scales = [scale / iteration for scale in standard_scales]\r\n","    return standard_scales, percss"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"1jB5PFsktvKe","executionInfo":{"status":"ok","timestamp":1612254056836,"user_tz":300,"elapsed":554448,"user":{"displayName":"Matthew So","photoUrl":"","userId":"18273672785190648211"}},"outputId":"65f4dcf6-0035-4cf7-c8b0-629f13360e96"},"source":["standard_scales, percss = nyul_train_dataloader(train_dataloader)"],"execution_count":null,"outputs":[{"output_type":"stream","text":["1\n","2\n","3\n","4\n","5\n","6\n","7\n","8\n","9\n","10\n","11\n","12\n","13\n","14\n","15\n","16\n","17\n","18\n","19\n","20\n","21\n","22\n","23\n","24\n","25\n","26\n","27\n","28\n","29\n","30\n","31\n","32\n","33\n","34\n","35\n","36\n","37\n","38\n","39\n","40\n","41\n","42\n","43\n","44\n","45\n","46\n","47\n","48\n","49\n","50\n","51\n","52\n","53\n","54\n","55\n","56\n","57\n","58\n","59\n","60\n","61\n","62\n","63\n","64\n","65\n","66\n","67\n","68\n","69\n","70\n","71\n","72\n","73\n","74\n","75\n","76\n","77\n","78\n","79\n","80\n","81\n","82\n","83\n","84\n","85\n","86\n","87\n","88\n","89\n","90\n","91\n","92\n","93\n","94\n","95\n","96\n","97\n","98\n","99\n","100\n","101\n","102\n","103\n","104\n","105\n","106\n","107\n","108\n","109\n","110\n","111\n","112\n","113\n","114\n","115\n","116\n","117\n","118\n","119\n","120\n","121\n","122\n","123\n","124\n","125\n","126\n","127\n","128\n","129\n","130\n","131\n","132\n","133\n","134\n","135\n","136\n","137\n","138\n","139\n","140\n","141\n","142\n","143\n","144\n","145\n","146\n","147\n","148\n","149\n","150\n","151\n","152\n","153\n","154\n","155\n","156\n","157\n","158\n","159\n","160\n","161\n","162\n","163\n","164\n","165\n","166\n","167\n","168\n","169\n","170\n","171\n","172\n","173\n","174\n","175\n","176\n","177\n","178\n","179\n","180\n","181\n","182\n","183\n","184\n","185\n","186\n","187\n","188\n","189\n","190\n","191\n","192\n","193\n","194\n","195\n","196\n","197\n","198\n","199\n","200\n","201\n","202\n","203\n","204\n","205\n","206\n","207\n","208\n","209\n","210\n","211\n","212\n","213\n","214\n","215\n","216\n","217\n","218\n","219\n","220\n","221\n","222\n","223\n","224\n","225\n","226\n","227\n","228\n","229\n","230\n","231\n","232\n","233\n","234\n","235\n","236\n","237\n","238\n","239\n","240\n","241\n","242\n","243\n","244\n","245\n","246\n","247\n","248\n","249\n","250\n","251\n","252\n","253\n","254\n","255\n","256\n","257\n","258\n","259\n","260\n","261\n","262\n","263\n","264\n","265\n","266\n","267\n","268\n","269\n","270\n","271\n","272\n","273\n","274\n","275\n","276\n","277\n","278\n","279\n","280\n","281\n","282\n","283\n","284\n","285\n","286\n","287\n","288\n","289\n","290\n","291\n","292\n","293\n","294\n"],"name":"stdout"}]},{"cell_type":"code","metadata":{"id":"pQ3tvV6MtZLY"},"source":["def dataloader_hist_norm(img_data, landmark_percs, standard_scale, seg_data):\r\n","    \"\"\"\r\n","    do the Nyul and Udupa histogram normalization routine with a given set of learned landmarks\r\n","    Args:\r\n","        img (nibabel.nifti1.Nifti1Image): image on which to find landmarks\r\n","        landmark_percs (np.ndarray): corresponding landmark points of standard scale\r\n","        standard_scale (np.ndarray): landmarks on the standard scale\r\n","        mask (nibabel.nifti1.Nifti1Image): foreground mask for img\r\n","    Returns:\r\n","        normalized (nibabel.nifti1.Nifti1Image): normalized image\r\n","    \"\"\"\r\n","    mask_data = seg_data\r\n","    mask_data[seg_data ==0] = 1\r\n","    mask_data = np.squeeze(mask_data, axis=0)\r\n","\r\n","    masked = img_data[mask_data > 0]\r\n","    landmarks = intensity_normalization.normalize.nyul.get_landmarks(masked, landmark_percs)\r\n","    f = interp1d(landmarks, standard_scale, fill_value='extrapolate')\r\n","    normed = f(img_data)\r\n","\r\n","    z = img_data\r\n","    z[img_data > 0] = normed[img_data > 0]\r\n","    return z #normed"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/","height":1000,"output_embedded_package_id":"1TnMQ4EcCGBnW64m9YxYGF3vyPxETitjQ"},"id":"yU6RWbGVsiv2","executionInfo":{"status":"error","timestamp":1612254198200,"user_tz":300,"elapsed":695804,"user":{"displayName":"Matthew So","photoUrl":"","userId":"18273672785190648211"}},"outputId":"8a424226-154b-48fb-ecc0-7d22a42107c7"},"source":["for all_img, seg_data in train_dataloader:\r\n","  for i, this_img in enumerate(all_img):\r\n","    if i == 0:\r\n","      transformed_img = dataloader_hist_norm(this_img, percss[i], standard_scales[i], seg_data)\r\n","      transformed_img = transformed_img[transformed_img>0]\r\n","      plt.hist(np.ravel(transformed_img), bins=30)\r\n","      plt.xlim(0, 150)\r\n","      plt.show()\r\n","     # plt.hist(np.ravel(this_img))\r\n","     # plt.show()"],"execution_count":null,"outputs":[{"output_type":"display_data","data":{"text/plain":"Output hidden; open in https://colab.research.google.com to view."},"metadata":{}}]}]}