{"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":"markdown","source":"This notebook uses the approach of generating bounding boxes from Class Activation Maps (CAM). CAM was introduced by Bolei Zhou et al. in the paper [Learning Deep Features for Discriminative Localization](https://arxiv.org/abs/1512.04150) and is a great tool for interpretation.","metadata":{}},{"cell_type":"code","source":"!pip install '../input/timm034/timm-0.3.4-py3-none-any.whl' -qq\n\n!cp /kaggle/input/gdcm-conda-install/gdcm.tar .\n!tar -xvzf gdcm.tar\n!conda install --offline ./gdcm/gdcm-2.8.9-py37h71b2a6d_0.tar.bz2\n!rm -rf ./gdcm.tar","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:57:19.674410Z","iopub.execute_input":"2021-07-20T00:57:19.675000Z","iopub.status.idle":"2021-07-20T00:58:23.367628Z","shell.execute_reply.started":"2021-07-20T00:57:19.674882Z","shell.execute_reply":"2021-07-20T00:58:23.366192Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import warnings\nwarnings.filterwarnings(\"ignore\", category=DeprecationWarning) ","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:58:36.561064Z","iopub.execute_input":"2021-07-20T00:58:36.561441Z","iopub.status.idle":"2021-07-20T00:58:36.567203Z","shell.execute_reply.started":"2021-07-20T00:58:36.561410Z","shell.execute_reply":"2021-07-20T00:58:36.565989Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from fastai.vision.all import *\nfrom fastai.medical.imaging import *\nimport pydicom\nimport cv2\nfrom timm import create_model\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut\nimport gdcm\nfrom torchvision.utils import save_image\nimport ast\nimport timm\nfrom glob import glob\nfrom skimage import exposure\n\nfrom fastai.vision.models.unet import _get_sz_change_idxs, hook_outputs\nmatplotlib.rcParams['image.cmap'] = 'viridis'","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:58:38.692552Z","iopub.execute_input":"2021-07-20T00:58:38.692916Z","iopub.status.idle":"2021-07-20T00:58:38.700714Z","shell.execute_reply.started":"2021-07-20T00:58:38.692887Z","shell.execute_reply":"2021-07-20T00:58:38.699318Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dicom_dataframe = pd.read_csv('../input/siim-covid/dicom_dataframe.csv')\n\nsource = '../input/siim-covid19-detection'\ndicom_dataframe = pd.read_csv('../input/siim-covid/dicom_dataframe.csv')\ntest_dataframe = pd.read_csv('../input/siim-covid/test_dataframe.csv')\nsub = pd.read_csv(f'{source}/sample_submission.csv')","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:58:40.673524Z","iopub.execute_input":"2021-07-20T00:58:40.673939Z","iopub.status.idle":"2021-07-20T00:58:41.197548Z","shell.execute_reply.started":"2021-07-20T00:58:40.673901Z","shell.execute_reply":"2021-07-20T00:58:41.196468Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Note in this case we will not be using the bounding box co-ordinates in the `train_image_level` csv file, we however have to extract the class labels for each image","metadata":{}},{"cell_type":"code","source":"train_image_level = pd.read_csv('../input/siim-covid19-detection/train_image_level.csv')\ntrain_study_level = pd.read_csv('../input/siim-covid19-detection/train_study_level.csv')\n\ntrain_study_level.rename(columns = {'id':'StudyInstanceUID'}, inplace = True)\n\ntrain_study_level['StudyInstanceUID'] = train_study_level['StudyInstanceUID'].apply(lambda x: f'{x[:12]}')\n\nmerged = pd.merge(train_image_level, train_study_level, on='StudyInstanceUID')\n\nmerged['id'] = merged['id'].apply(lambda x: f'{x[:-6]}')\nmerged.loc[merged['Negative for Pneumonia']==1, 'class_y'] = 'negative'\nmerged.loc[merged['Typical Appearance']==1, 'class_y'] = 'typical'\nmerged.loc[merged['Indeterminate Appearance']==1, 'class_y'] = 'indeterminate'\nmerged.loc[merged['Atypical Appearance']==1, 'class_y'] = 'atypical'\nmerged.drop(['boxes', 'Negative for Pneumonia', 'Typical Appearance', \n             'Indeterminate Appearance', 'Atypical Appearance'], axis=1, inplace=True)\ndel merged['label']","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:58:46.412334Z","iopub.execute_input":"2021-07-20T00:58:46.412698Z","iopub.status.idle":"2021-07-20T00:58:46.532752Z","shell.execute_reply.started":"2021-07-20T00:58:46.412669Z","shell.execute_reply":"2021-07-20T00:58:46.531547Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"merged.rename(columns = {'id':'SOPInstanceUID'}, inplace = True)\nmerged[:2]","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:58:47.339577Z","iopub.execute_input":"2021-07-20T00:58:47.339937Z","iopub.status.idle":"2021-07-20T00:58:47.359649Z","shell.execute_reply.started":"2021-07-20T00:58:47.339906Z","shell.execute_reply":"2021-07-20T00:58:47.358840Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Merge with the dataframe which contains all the dicom metadata.  You can see this [notebook](https://www.kaggle.com/avirdee/siim-covid-19-initial-pipeline-fastai) for learn more about creating the dataframe.  We now have the file paths to each image which we will need when creating the datablock.","metadata":{}},{"cell_type":"code","source":"dicom_merge = pd.merge(dicom_dataframe, merged, on='SOPInstanceUID')","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:58:48.421618Z","iopub.execute_input":"2021-07-20T00:58:48.422182Z","iopub.status.idle":"2021-07-20T00:58:48.459540Z","shell.execute_reply.started":"2021-07-20T00:58:48.422119Z","shell.execute_reply":"2021-07-20T00:58:48.458482Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = dicom_merge[['SOPInstanceUID', 'class_y', 'fname']].copy()\ndf[:5]","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:58:48.795273Z","iopub.execute_input":"2021-07-20T00:58:48.795666Z","iopub.status.idle":"2021-07-20T00:58:48.817883Z","shell.execute_reply.started":"2021-07-20T00:58:48.795633Z","shell.execute_reply":"2021-07-20T00:58:48.816678Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":" create a new class` FixMonochrome` that inherits from `PILDicom` to fix the monochrome issue.","metadata":{}},{"cell_type":"code","source":"class FixMonochrome(PILDicom):\n    @classmethod\n    def create(cls, fn:(Path, str, bytes), fix_monochrome = True)->None:\n        if isinstance(fn, bytes): im = pydicom.dcmread(pydicom.filebase.DicomBytesIO(fn))\n        if isinstance(fn, (Path, str)): im = pydicom.dcmread(fn)\n        scaled = np.array(im.pixel_array)   \n        if fix_monochrome and im.PhotometricInterpretation == \"MONOCHROME1\":\n            scaled = np.amax(scaled) - scaled\n        scaled = scaled - np.min(scaled)\n        scaled = scaled / np.max(scaled)\n        scaled = (scaled * 255).astype(np.uint8)\n        pill_im = Image.fromarray(scaled)\n        return cls(pill_im)","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:58:49.502033Z","iopub.execute_input":"2021-07-20T00:58:49.502421Z","iopub.status.idle":"2021-07-20T00:58:49.510358Z","shell.execute_reply.started":"2021-07-20T00:58:49.502390Z","shell.execute_reply":"2021-07-20T00:58:49.509549Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"batch_tfms = [Resize(128), *aug_transforms(do_flip=False, \n                                           flip_vert=False, \n                                           max_rotate=17,  \n                                           min_zoom=1.,\n                                           max_zoom=1.7,\n                                           max_lighting=0.1), Normalize.from_stats(*imagenet_stats)]","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:58:50.103047Z","iopub.execute_input":"2021-07-20T00:58:50.103593Z","iopub.status.idle":"2021-07-20T00:58:50.115109Z","shell.execute_reply.started":"2021-07-20T00:58:50.103559Z","shell.execute_reply":"2021-07-20T00:58:50.114129Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"siim = DataBlock(blocks=(ImageBlock(cls=FixMonochrome), CategoryBlock),\n                 get_x=ColReader('fname'),\n                 splitter=RandomSplitter(),\n                 get_y=ColReader('class_y'),\n                 item_tfms=Resize(264))\n\ndls = siim.dataloaders(df, bs=32)\ndls.show_batch(max_n=16, nrows=2)","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:58:51.839225Z","iopub.execute_input":"2021-07-20T00:58:51.839606Z","iopub.status.idle":"2021-07-20T00:59:10.752889Z","shell.execute_reply.started":"2021-07-20T00:58:51.839572Z","shell.execute_reply":"2021-07-20T00:59:10.751975Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Using `resnet200d`","metadata":{}},{"cell_type":"code","source":"os.makedirs('/root/.cache/torch/hub/checkpoints/')\n!cp '../input/timm034/resnet200d_ra2-bdba9bf9.pth' '/root/.cache/torch/hub/checkpoints/resnet200d_ra2-bdba9bf9.pth'","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:59:10.754065Z","iopub.execute_input":"2021-07-20T00:59:10.754531Z","iopub.status.idle":"2021-07-20T00:59:14.306981Z","shell.execute_reply.started":"2021-07-20T00:59:10.754499Z","shell.execute_reply":"2021-07-20T00:59:14.305662Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To use [timm](https://github.com/rwightman/pytorch-image-models) models with fastai we use the code from this [notebook](https://github.com/muellerzr/Practical-Deep-Learning-for-Coders-2.0/blob/master/Computer%20Vision/05_EfficientNet_and_Custom_Weights.ipynb)","metadata":{}},{"cell_type":"code","source":"def create_timm_body(arch:str, pretrained=True, cut=None):\n    model = create_model(arch, pretrained=pretrained)\n    if cut is None:\n        ll = list(enumerate(model.children()))\n        cut = next(i for i,o in reversed(ll) if has_pool_type(o))\n    if isinstance(cut, int): return nn.Sequential(*list(model.children())[:cut])\n    elif callable(cut): return cut(model)\n    else: raise NamedError(\"cut must be either integer or function\")\n        \nbody = create_timm_body('resnet200d', pretrained=True)\n\nl = nn.Conv2d(1, 64, kernel_size=(7,7), stride=(2,2),\n                    padding=(3,3), bias=False)\nl.weight = nn.Parameter(l.weight.sum(dim=1, keepdim=True))\n\nbody[0] = l\n\nhead = create_head(2048, dls.c)\nmodel = nn.Sequential(body, head)\napply_init(model[1], nn.init.kaiming_normal_)","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:59:14.309264Z","iopub.execute_input":"2021-07-20T00:59:14.309599Z","iopub.status.idle":"2021-07-20T00:59:16.128967Z","shell.execute_reply.started":"2021-07-20T00:59:14.309567Z","shell.execute_reply":"2021-07-20T00:59:16.127891Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Create the `Learner`","metadata":{}},{"cell_type":"code","source":"learn = Learner(dls, \n                model, \n                loss_func=LabelSmoothingCrossEntropyFlat(),\n                metrics=accuracy,\n                cbs=[ShowGraphCallback()])","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:59:16.130502Z","iopub.execute_input":"2021-07-20T00:59:16.130775Z","iopub.status.idle":"2021-07-20T00:59:16.137575Z","shell.execute_reply.started":"2021-07-20T00:59:16.130748Z","shell.execute_reply":"2021-07-20T00:59:16.136353Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Load an earlier trained ` model`","metadata":{}},{"cell_type":"code","source":"Path('./models').mkdir(exist_ok=True, parents=True)\n\n!cp '../input/siim-covid/siim_two.pth' './models/siim_two.pth'","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:59:16.138897Z","iopub.execute_input":"2021-07-20T00:59:16.139259Z","iopub.status.idle":"2021-07-20T00:59:24.265021Z","shell.execute_reply.started":"2021-07-20T00:59:16.139211Z","shell.execute_reply":"2021-07-20T00:59:24.263825Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"learn.load('siim_two')","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:59:24.266991Z","iopub.execute_input":"2021-07-20T00:59:24.267318Z","iopub.status.idle":"2021-07-20T00:59:25.221401Z","shell.execute_reply.started":"2021-07-20T00:59:24.267286Z","shell.execute_reply":"2021-07-20T00:59:25.220377Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Create a new dataframe and add the file paths of the images in the test dataset to help with inference","metadata":{}},{"cell_type":"code","source":"study_level = sub[sub['id'].str.endswith('_study')]\nimage_level = sub[sub['id'].str.endswith('_image')]\n\nimage_level['id2'] = image_level['id']\nimage_level['id2'] = image_level['id2'].apply(lambda x: f'{x[:-6]}')\nimage_level.rename(columns = {'id2':'SOPInstanceUID'}, inplace = True)\ntest_image_merge = pd.merge(test_dataframe, image_level, on='SOPInstanceUID')\ntest_img = test_image_merge[['id', 'fname']].copy()\ntest_img[:5]","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:59:25.222832Z","iopub.execute_input":"2021-07-20T00:59:25.223259Z","iopub.status.idle":"2021-07-20T00:59:25.268781Z","shell.execute_reply.started":"2021-07-20T00:59:25.223215Z","shell.execute_reply":"2021-07-20T00:59:25.267752Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Grad CAM","metadata":{}},{"cell_type":"markdown","source":"**CAM** uses the output of the last convolutional layer to provide a heatmap visualization.\n\nTo do this we need a way to access the activations during the training phase and this in `Pytorch` is done using a `hook`. A `hook` is basically a function that is executed when either the `forward` or `backward` pass is called.\n\n`Fastai` has a handy `Hook` class that accesses the activations in the forward pass","metadata":{}},{"cell_type":"code","source":"class Hook():\n    def __init__(self, m):\n        self.hook = m.register_forward_hook(self.hook_func)   \n    def hook_func(self, m, i, o): self.stored = o.detach().clone()\n    def __enter__(self, *args): return self\n    def __exit__(self, *args): self.hook.remove()","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:59:25.271841Z","iopub.execute_input":"2021-07-20T00:59:25.272485Z","iopub.status.idle":"2021-07-20T00:59:25.279861Z","shell.execute_reply.started":"2021-07-20T00:59:25.272425Z","shell.execute_reply":"2021-07-20T00:59:25.278567Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Grad CAM** is a variation of CAM that was introduced in [Grad-CAM: Why Did You Say That? Visual Explanations from Deep Networks via Gradient-based Localization](https://arxiv.org/abs/1611.07450).\n\nGrad CAM is similar to CAM in that it is a weighted combination of feature maps but now followed by a `ReLU`. The output of Grad CAM is a “class-discriminative localization map” where the hot part of the heatmap corresponds to a particular class.  The gradients of every layer are calculated by PyTorch during the backward pass.  To access the gradients we can register a hook on the backward pass and store these gradients.  `Fastai` has a `HookBwd` class that can intercept these stored gradients from the backward pass.","metadata":{}},{"cell_type":"code","source":"class HookBwd():\n    def __init__(self, m):\n        self.hook = m.register_backward_hook(self.hook_func)   \n    def hook_func(self, m, gi, go): self.stored = go[0].detach().clone()\n    def __enter__(self, *args): return self\n    def __exit__(self, *args): self.hook.remove()","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:59:25.281728Z","iopub.execute_input":"2021-07-20T00:59:25.282102Z","iopub.status.idle":"2021-07-20T00:59:25.294243Z","shell.execute_reply.started":"2021-07-20T00:59:25.282066Z","shell.execute_reply":"2021-07-20T00:59:25.293330Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The great thing about Grad CAM is that we will be able to access the heatmaps for any of the classes in this dataset.\n\nIn this case there are 4 classes","metadata":{}},{"cell_type":"code","source":"dls.vocab","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:59:25.295467Z","iopub.execute_input":"2021-07-20T00:59:25.295943Z","iopub.status.idle":"2021-07-20T00:59:25.310939Z","shell.execute_reply.started":"2021-07-20T00:59:25.295883Z","shell.execute_reply":"2021-07-20T00:59:25.309970Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### How do we generate heatmaps?","metadata":{}},{"cell_type":"markdown","source":"Grab a test file from the test dataframe created earlier","metadata":{}},{"cell_type":"code","source":"im = test_img[300:301]['fname'].item()","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:59:25.312534Z","iopub.execute_input":"2021-07-20T00:59:25.313134Z","iopub.status.idle":"2021-07-20T00:59:25.322945Z","shell.execute_reply.started":"2021-07-20T00:59:25.313095Z","shell.execute_reply":"2021-07-20T00:59:25.321927Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can use `learn.predict` to predict the class of the test image.  This returns the `class`, `prediction` and `probabilities`","metadata":{}},{"cell_type":"code","source":"cl, pred, probs = learn.predict(im)\nprint(cl, pred, probs)","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:59:25.324459Z","iopub.execute_input":"2021-07-20T00:59:25.325054Z","iopub.status.idle":"2021-07-20T00:59:26.551447Z","shell.execute_reply.started":"2021-07-20T00:59:25.325019Z","shell.execute_reply":"2021-07-20T00:59:26.550610Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To be able to generate the correct heatmap for the dominant class `typical` in this case, we need to pass the `prediction` value which is `3` when generating the heatmap:","metadata":{}},{"cell_type":"code","source":"cls = pred.item()","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:59:26.554200Z","iopub.execute_input":"2021-07-20T00:59:26.556101Z","iopub.status.idle":"2021-07-20T00:59:26.563816Z","shell.execute_reply.started":"2021-07-20T00:59:26.556052Z","shell.execute_reply":"2021-07-20T00:59:26.562659Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"x, = first(dls.test_dl([im]))\nx_dec = TensorImage(dls.train.decode((x,))[0][0])","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:59:26.570747Z","iopub.execute_input":"2021-07-20T00:59:26.571111Z","iopub.status.idle":"2021-07-20T00:59:27.011356Z","shell.execute_reply.started":"2021-07-20T00:59:26.571080Z","shell.execute_reply":"2021-07-20T00:59:27.010061Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Grab the gradients using `Hook` and `HookBwd` and here we are grabbing them from the last layer.","metadata":{}},{"cell_type":"code","source":"with HookBwd(learn.model[0][-1]) as hookg:\n    with Hook(learn.model[0][-1]) as hook:\n        output = learn.model.eval()(x.cpu())\n        act = hook.stored\n    output[0, cls].backward()\n    grad = hookg.stored","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:59:27.012993Z","iopub.execute_input":"2021-07-20T00:59:27.013415Z","iopub.status.idle":"2021-07-20T00:59:29.019610Z","shell.execute_reply.started":"2021-07-20T00:59:27.013374Z","shell.execute_reply":"2021-07-20T00:59:29.018423Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Generate the heatmap","metadata":{}},{"cell_type":"code","source":"w = grad[0].mean(dim=[1,2], keepdim=True)\ncam_map = (w * act[0]).sum(0)","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:59:29.021194Z","iopub.execute_input":"2021-07-20T00:59:29.021542Z","iopub.status.idle":"2021-07-20T00:59:29.029430Z","shell.execute_reply.started":"2021-07-20T00:59:29.021509Z","shell.execute_reply":"2021-07-20T00:59:29.028107Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Check the shape of `cam_map`","metadata":{}},{"cell_type":"code","source":"cam_map.shape","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:59:29.030829Z","iopub.execute_input":"2021-07-20T00:59:29.031252Z","iopub.status.idle":"2021-07-20T00:59:29.045639Z","shell.execute_reply.started":"2021-07-20T00:59:29.031195Z","shell.execute_reply":"2021-07-20T00:59:29.044699Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In this case the heatmap is a 9 by 9 square.  Note that the shape of `cam_map` is determined by the `item_tfms` shape set in the `DataBlock`.  The larger the image used during training the larger the heat map.\n\nWe can view the heatmap using `show_image`","metadata":{}},{"cell_type":"code","source":"show_image(cam_map, figsize=(7,7));","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:59:29.046845Z","iopub.execute_input":"2021-07-20T00:59:29.047265Z","iopub.status.idle":"2021-07-20T00:59:29.164342Z","shell.execute_reply.started":"2021-07-20T00:59:29.047232Z","shell.execute_reply":"2021-07-20T00:59:29.163191Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The location on the heatmap that has the highest activations is the bright yellow square at position (3, 3).\n\nWe can view the heatmap on the image like so:","metadata":{}},{"cell_type":"code","source":"_,ax = plt.subplots()\nx_dec.show(ctx=ax)\nax.imshow(cam_map.detach().cpu(), alpha=0.6, extent=(0, x_dec.shape[1], x_dec.shape[2],0),\n              interpolation='bilinear', cmap='magma');","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:59:41.545256Z","iopub.execute_input":"2021-07-20T00:59:41.545674Z","iopub.status.idle":"2021-07-20T00:59:41.683576Z","shell.execute_reply.started":"2021-07-20T00:59:41.545639Z","shell.execute_reply":"2021-07-20T00:59:41.682572Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can see that the location with the highest activations is just left of the sternum bone.  There are 2 other locations that have activations but I am currently only interested in the location with the highest activations.","metadata":{}},{"cell_type":"markdown","source":"Now we need to be able to generate a bounding box around the area with the most activations.\n\nTo do this we need to iterate over the `cam_map` values and get the indexes of the x and and y locations where the activations are the highest.\n\n`cam_map` is a 9 by 9 `TensorDicom`","metadata":{}},{"cell_type":"code","source":"cam_map, cam_map.shape","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:59:49.723260Z","iopub.execute_input":"2021-07-20T00:59:49.723806Z","iopub.status.idle":"2021-07-20T00:59:49.731499Z","shell.execute_reply.started":"2021-07-20T00:59:49.723771Z","shell.execute_reply":"2021-07-20T00:59:49.730533Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cms = cam_map.shape[0]","metadata":{"execution":{"iopub.status.busy":"2021-07-20T00:59:52.784101Z","iopub.execute_input":"2021-07-20T00:59:52.784531Z","iopub.status.idle":"2021-07-20T00:59:52.788476Z","shell.execute_reply.started":"2021-07-20T00:59:52.784495Z","shell.execute_reply":"2021-07-20T00:59:52.787651Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Convert it to a tensor","metadata":{}},{"cell_type":"code","source":"t = np.array(cam_map)\ntr = tensor(t)","metadata":{"execution":{"iopub.status.busy":"2021-07-20T01:00:03.376862Z","iopub.execute_input":"2021-07-20T01:00:03.377437Z","iopub.status.idle":"2021-07-20T01:00:03.381938Z","shell.execute_reply.started":"2021-07-20T01:00:03.377400Z","shell.execute_reply":"2021-07-20T01:00:03.381216Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"What is the maximum value in the whole tensor?","metadata":{}},{"cell_type":"code","source":"tr.max()","metadata":{"execution":{"iopub.status.busy":"2021-07-20T01:00:09.632245Z","iopub.execute_input":"2021-07-20T01:00:09.632605Z","iopub.status.idle":"2021-07-20T01:00:09.640111Z","shell.execute_reply.started":"2021-07-20T01:00:09.632577Z","shell.execute_reply":"2021-07-20T01:00:09.639040Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To get the indexes of the maximum tensor we start by first finding the maximum value in each row","metadata":{}},{"cell_type":"code","source":"val = []\nfor i in range(0, tr.shape[0]):\n    index, value = max(enumerate(tr[i]), key=operator.itemgetter(1))\n    val.append(value)","metadata":{"execution":{"iopub.status.busy":"2021-07-20T01:08:02.959547Z","iopub.execute_input":"2021-07-20T01:08:02.959980Z","iopub.status.idle":"2021-07-20T01:08:02.967429Z","shell.execute_reply.started":"2021-07-20T01:08:02.959942Z","shell.execute_reply":"2021-07-20T01:08:02.966572Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"val","metadata":{"execution":{"iopub.status.busy":"2021-07-20T01:08:13.890682Z","iopub.execute_input":"2021-07-20T01:08:13.891053Z","iopub.status.idle":"2021-07-20T01:08:13.901046Z","shell.execute_reply.started":"2021-07-20T01:08:13.891022Z","shell.execute_reply":"2021-07-20T01:08:13.899906Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We now have the maximum value of each row.  Now to get the y_index and confirm that it matches the maximum value","metadata":{}},{"cell_type":"code","source":"y_index, y_value = max(enumerate(val), key=operator.itemgetter(1))\nprint(y_index, y_value)","metadata":{"execution":{"iopub.status.busy":"2021-07-20T01:10:09.852045Z","iopub.execute_input":"2021-07-20T01:10:09.852486Z","iopub.status.idle":"2021-07-20T01:10:09.859316Z","shell.execute_reply.started":"2021-07-20T01:10:09.852451Z","shell.execute_reply":"2021-07-20T01:10:09.858476Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Get the x_index","metadata":{}},{"cell_type":"code","source":"x_index, x_value = max(enumerate(tr[y_index]), key=operator.itemgetter(1))\nx_index, x_value","metadata":{"execution":{"iopub.status.busy":"2021-07-20T01:10:19.983806Z","iopub.execute_input":"2021-07-20T01:10:19.984363Z","iopub.status.idle":"2021-07-20T01:10:19.991111Z","shell.execute_reply.started":"2021-07-20T01:10:19.984328Z","shell.execute_reply":"2021-07-20T01:10:19.990335Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The maximum value of `0.0261` can be found at x_index 2 and y_index 5\n\nNow we need to match that index to the image\n\nTo do this we need to know the shape of the image","metadata":{}},{"cell_type":"code","source":"imz = pydicom.dcmread(im).pixel_array \nimz.shape","metadata":{"execution":{"iopub.status.busy":"2021-07-20T01:10:54.353818Z","iopub.execute_input":"2021-07-20T01:10:54.354194Z","iopub.status.idle":"2021-07-20T01:10:54.401590Z","shell.execute_reply.started":"2021-07-20T01:10:54.354160Z","shell.execute_reply":"2021-07-20T01:10:54.400396Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"And since the heatmap is broken down into 9 by 9 squares we need to divide the image into 9 by 9 sections so that they correspond to each section of the heatmap.\n\nAlong the x-axis","metadata":{}},{"cell_type":"code","source":"x_ = imz.shape[0] // cms\nx_","metadata":{"execution":{"iopub.status.busy":"2021-07-20T01:11:51.496006Z","iopub.execute_input":"2021-07-20T01:11:51.496647Z","iopub.status.idle":"2021-07-20T01:11:51.502362Z","shell.execute_reply.started":"2021-07-20T01:11:51.496608Z","shell.execute_reply":"2021-07-20T01:11:51.501554Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Along the y-axis","metadata":{}},{"cell_type":"code","source":"y_ = imz.shape[1] // cms\ny_","metadata":{"execution":{"iopub.status.busy":"2021-07-20T01:11:54.768383Z","iopub.execute_input":"2021-07-20T01:11:54.768890Z","iopub.status.idle":"2021-07-20T01:11:54.774608Z","shell.execute_reply.started":"2021-07-20T01:11:54.768857Z","shell.execute_reply":"2021-07-20T01:11:54.773733Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So the image now has 9 by 9 squares each of size 332 by 332.\n\nTo match the square with highest activations we multiply the indexes we calculated earlier to the size of each square","metadata":{}},{"cell_type":"code","source":"xx = x_index * x_\nxx","metadata":{"execution":{"iopub.status.busy":"2021-07-20T01:11:59.056273Z","iopub.execute_input":"2021-07-20T01:11:59.056666Z","iopub.status.idle":"2021-07-20T01:11:59.063240Z","shell.execute_reply.started":"2021-07-20T01:11:59.056635Z","shell.execute_reply":"2021-07-20T01:11:59.061612Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"yy = y_index * y_\nyy","metadata":{"execution":{"iopub.status.busy":"2021-07-20T01:12:05.079968Z","iopub.execute_input":"2021-07-20T01:12:05.080349Z","iopub.status.idle":"2021-07-20T01:12:05.086746Z","shell.execute_reply.started":"2021-07-20T01:12:05.080318Z","shell.execute_reply":"2021-07-20T01:12:05.085549Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Using `TensorBBox` we can now plot a bounding box around the area of the image which has the highest activations","metadata":{}},{"cell_type":"code","source":"box = TensorBBox([xx,yy, 600 + xx, yy - 600])\n_,ax = plt.subplots()\nx_dec.show(ctx=ax)\nax.imshow(cam_map.detach().cpu(), alpha=0.6, extent=(0, x_dec.shape[1], x_dec.shape[2],0),\n              interpolation='bilinear', cmap='magma');\nctx = show_image(imz)\nbox.show(ctx=ctx);","metadata":{"execution":{"iopub.status.busy":"2021-07-20T01:13:32.131631Z","iopub.execute_input":"2021-07-20T01:13:32.132078Z","iopub.status.idle":"2021-07-20T01:13:34.403864Z","shell.execute_reply.started":"2021-07-20T01:13:32.132040Z","shell.execute_reply":"2021-07-20T01:13:34.403055Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_b(im):\n    cl, pred, probs = learn.predict(im)\n    cls = pred.item()\n    x, = first(dls.test_dl([im]))\n    x_dec = TensorImage(dls.train.decode((x,))[0][0])\n    with HookBwd(learn.model[0][-1]) as hookg:\n        with Hook(learn.model[0][-1]) as hook:\n            output = learn.model.eval()(x.cpu())\n            act = hook.stored\n        output[0, cls].backward()\n        grad = hookg.stored\n    w = grad[0].mean(dim=[1,2], keepdim=True)\n    cam_map = (w * act[0]).sum(0)\n    t = np.array(cam_map)\n    tr = tensor(t)\n    val = []\n    for i in range(0, tr.shape[0]):\n        index, value = max(enumerate(tr[i]), key=operator.itemgetter(1))\n        val.append(value)\n        \n    y_index, y_value = max(enumerate(val), key=operator.itemgetter(1))\n    x_index, x_value = max(enumerate(tr[y_index]), key=operator.itemgetter(1))\n    \n    imz = pydicom.dcmread(im).pixel_array \n    x_ = imz.shape[0] // cms\n    y_ = imz.shape[1] // cms\n    xx = x_index * x_\n    yy = y_index * y_\n    \n    box = TensorBBox([xx,yy, x_ + xx, y_ + yy])\n    _,ax = plt.subplots()\n    x_dec.show(ctx=ax)\n    ax.imshow(cam_map.detach().cpu(), alpha=0.6, extent=(0, x_dec.shape[1], x_dec.shape[2],0),\n              interpolation='bilinear', cmap='magma');\n    ctx = show_image(imz)\n    box.show(ctx=ctx);","metadata":{"execution":{"iopub.status.busy":"2021-07-20T01:34:43.317530Z","iopub.execute_input":"2021-07-20T01:34:43.317984Z","iopub.status.idle":"2021-07-20T01:34:43.723592Z","shell.execute_reply.started":"2021-07-20T01:34:43.317942Z","shell.execute_reply":"2021-07-20T01:34:43.722280Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"im = test_img[201:202]['fname'].item()\ncl, pred, probs = learn.predict(im)\ncls = pred.item()\nx, = first(dls.test_dl([im]))\nx_dec = TensorImage(dls.train.decode((x,))[0][0])\nwith HookBwd(learn.model[0][-1]) as hookg:\n    with Hook(learn.model[0][-1]) as hook:\n        output = learn.model.eval()(x.cpu())\n        act = hook.stored\n    output[0, cls].backward()\n    grad = hookg.stored\nw = grad[0].mean(dim=[1,2], keepdim=True)\ncam_map = (w * act[0]).sum(0)\nt = np.array(cam_map)\ntr = tensor(t)\nval = []\nfor i in range(0, tr.shape[0]):\n    index, value = max(enumerate(tr[i]), key=operator.itemgetter(1))\n    val.append(value)\n        \ny_index, y_value = max(enumerate(val), key=operator.itemgetter(1))\nx_index, x_value = max(enumerate(tr[y_index]), key=operator.itemgetter(1))\n\nprint (val)\nprint(x_index, y_index)","metadata":{"execution":{"iopub.status.busy":"2021-07-20T01:36:58.692994Z","iopub.execute_input":"2021-07-20T01:36:58.693420Z","iopub.status.idle":"2021-07-20T01:37:01.509737Z","shell.execute_reply.started":"2021-07-20T01:36:58.693373Z","shell.execute_reply":"2021-07-20T01:37:01.508404Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"imz = pydicom.dcmread(im).pixel_array \nprint(imz.shape)","metadata":{"execution":{"iopub.status.busy":"2021-07-20T01:55:19.235352Z","iopub.execute_input":"2021-07-20T01:55:19.235915Z","iopub.status.idle":"2021-07-20T01:55:19.258512Z","shell.execute_reply.started":"2021-07-20T01:55:19.235871Z","shell.execute_reply":"2021-07-20T01:55:19.256868Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"imz = pydicom.dcmread(im).pixel_array \nprint(imz.shape)\nx_ = imz.shape[1] // cms\nprint(x_)\ny_ = imz.shape[1] // cms\nprint(y_)\nxx = x_index * x_\nyy = y_index * y_","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"xx = x_index * x_\nyy = y_index * y_\nprint(xx, yy)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"box = TensorBBox([xx,yy, xx+10,yy+10])\n_,ax = plt.subplots()\nx_dec.show(ctx=ax)\nax.imshow(cam_map.detach().cpu(), alpha=0.6, extent=(0, x_dec.shape[1], x_dec.shape[2],0),\n            interpolation='bilinear', cmap='magma');\nctx = show_image(imz)\nbox.show(ctx=ctx);","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i in range(201,202):\n    testi = test_img[i:i+1]['fname'].item()\n    get_b(testi)","metadata":{"execution":{"iopub.status.busy":"2021-07-20T01:34:46.358696Z","iopub.execute_input":"2021-07-20T01:34:46.359101Z","iopub.status.idle":"2021-07-20T01:34:49.876330Z","shell.execute_reply.started":"2021-07-20T01:34:46.359069Z","shell.execute_reply":"2021-07-20T01:34:49.874636Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We now have a bounding box around the area with the highest activations","metadata":{}}]}