{"cells":[{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["# Imports\n","import tensorflow as tf\n","import cv2\n","import numpy as np\n","import pandas as pd\n","import matplotlib.pyplot as plt\n","from tensorflow.keras.preprocessing.image import load_img\n","from PIL import Image\n","from scipy import ndimage\n","import matplotlib.patches as patches\n","from skimage.filters import threshold_multiotsu\n","from cv2 import imshow as cv2_imshow\n","import os\n","import matplotlib.patches as patches"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["# To visualize which parts of the image are the most important for the each \n","#label for the classifier, let’s set up the Grad-CAM process\n","def xray_CAM(model,img_array, last_conv_layer_name, pred_index): \n"," # First, we create a model that maps the input image to the activations\n"," # of the last conv layer as well as the output predictions\n"," grad_model = tf.keras.models.Model(\n"," [model.inputs], [model.get_layer(last_conv_layer_name).output, model.output]\n"," )\n","\n"," # Then, we compute the gradient of the top predicted class for our input image\n"," # with respect to the activations of the last conv layer\n"," with tf.GradientTape() as tape:\n"," last_conv_layer_output, preds = grad_model(img_array)\n"," class_channel = preds[:, pred_index]\n","\n"," # This is the gradient of the output neuron (top predicted or chosen)\n"," # with regard to the output feature map of the last conv layer\n"," grads = tape.gradient(class_channel, last_conv_layer_output)\n","\n"," # This is a vector where each entry is the mean intensity of the gradient\n"," # over a specific feature map channel\n"," pooled_grads = tf.reduce_mean(grads, axis=(0, 1, 2))\n","\n"," # We multiply each channel in the feature map array\n"," # by \"how important this channel is\" with regard to the top predicted class\n"," # then sum all the channels to obtain the heatmap class activation\n"," last_conv_layer_output = last_conv_layer_output[0]\n"," heatmap = last_conv_layer_output @ pooled_grads[..., tf.newaxis]\n"," heatmap = tf.squeeze(heatmap)\n","\n"," # For visualization purpose, we will also normalize the heatmap between 0 & 1\n"," heatmap = tf.maximum(heatmap, 0) / tf.math.reduce_max(heatmap)\n"," return heatmap.numpy()\n","\n","\n","def superimpose(img_path, heatmap):\n"," # We use cv2 to load the original image\n"," img = cv2.imread(img_path)\n","\n"," # We resize the heatmap to have the same size as the original image\n"," heatmap = cv2.resize(heatmap, (img.shape[1], img.shape[0]))\n","\n"," # We convert the heatmap to RGB\n"," heatmap = np.uint8(255 * heatmap)\n","\n"," # We apply the heatmap to the original image\n"," heatmap = cv2.applyColorMap(heatmap, cv2.COLORMAP_JET)\n","\n"," # 0.4 here is a heatmap intensity factor\n"," superimposed_img = heatmap * 0.4 + img\n","\n"," ## Save the image to disk\n"," # results folder must exist\n"," cv2.imwrite('results/xray_cam.jpg', superimposed_img)\n"," img_samp = cv2.imread('results/xray_cam.jpg') #OpenCV has BGR order\n"," img_samp = cv2.cvtColor(img_samp, cv2.COLOR_BGR2RGB) #matplotlib has RGB order\n"," \n"," return img_samp\n"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["def grad_cam_plus(model, img,\n"," layer_name,\n"," category_id):\n"," conv_layer = model.get_layer(layer_name)\n"," heatmap_model = tf.keras.models.Model([model.inputs], [conv_layer.output, model.output])\n","\n"," with tf.GradientTape() as gtape1:\n"," with tf.GradientTape() as gtape2:\n"," with tf.GradientTape() as gtape3:\n"," conv_output, predictions = heatmap_model(img)\n"," if category_id==None:\n"," category_id = np.argmax(predictions[0])\n"," output = predictions[:, category_id]\n"," conv_first_grad = gtape3.gradient(output, conv_output)\n"," conv_second_grad = gtape2.gradient(conv_first_grad, conv_output)\n"," conv_third_grad = gtape1.gradient(conv_second_grad, conv_output)\n","\n"," global_sum = np.sum(conv_output, axis=(0, 1, 2))\n","\n"," alpha_num = conv_second_grad[0]\n"," alpha_denom = conv_second_grad[0]*2.0 + conv_third_grad[0]*global_sum\n"," alpha_denom = np.where(alpha_denom != 0.0, alpha_denom, 1e-10)\n","\n"," \n"," alphas = alpha_num/alpha_denom\n"," alpha_normalization_constant = np.sum(alphas, axis=(0,1))\n"," alphas /= alpha_normalization_constant\n","\n","\n"," weights = np.maximum(conv_first_grad[0], 0.0)\n","\n","\n"," deep_linearization_weights = np.sum(weights*alphas, axis=(0,1))\n"," deep_linearization_weights[np.isnan(deep_linearization_weights)] = -1\n","\n"," cam_map=deep_linearization_weights*conv_output[0]\n","\n"," grad_CAM_map = np.sum(cam_map, axis=2)\n","\n","\n","\n"," heatmap = np.maximum(grad_CAM_map, 0)\n"," max_heat = np.max(heatmap)\n"," if max_heat == 0:\n"," max_heat = 1e-10\n"," heatmap /= max_heat\n","\n"," return heatmap"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["def ScoreCam(model, img_array, layer_name, cls, max_N=-1, si=True):\n","\n"," if cls==-1:\n"," cls = np.argmax(model.predict(img_array))\n"," act_map_array = tf.keras.models.Model(inputs=model.input, outputs=model.get_layer(layer_name).output).predict(img_array)\n"," \n"," # extract effective maps\n"," if max_N != -1:\n"," act_map_std_list = [np.std(act_map_array[0,:,:,k]) for k in range(act_map_array.shape[3])]\n"," unsorted_max_indices = np.argpartition(-np.array(act_map_std_list), max_N)[:max_N]\n"," max_N_indices = unsorted_max_indices[np.argsort(-np.array(act_map_std_list)[unsorted_max_indices])]\n"," act_map_array = act_map_array[:,:,:,max_N_indices]\n","\n"," input_shape = model.layers[0].output_shape[0][1:] # get input shape\n"," # 1. upsample to original input size\n"," act_map_resized_list = [np.resize(act_map_array[0,:,:,k], input_shape[:2]) for k in range(act_map_array.shape[3])]\n"," # 2. normalize the raw activation value in each activation map into [0, 1]\n"," act_map_normalized_list = []\n"," for act_map_resized in act_map_resized_list:\n"," if np.max(act_map_resized) - np.min(act_map_resized) != 0:\n"," if si:\n"," act_map_normalized = (act_map_resized - np.min(act_map_resized)) / (np.max(act_map_resized) - np.min(act_map_resized))\n"," else:\n"," act_map_normalized = (act_map_resized - np.mean(act_map_resized)) / (np.max(act_map_resized) - np.min(act_map_resized))\n"," else:\n"," act_map_normalized = act_map_resized\n"," act_map_normalized_list.append(act_map_normalized)\n"," # 3. project highlighted area in the activation map to original input space by multiplying the normalized activation map\n"," masked_input_list = []\n"," #print(img_array)\n"," for act_map_normalized in act_map_normalized_list:\n"," masked_input = np.copy(img_array)\n"," for k in range(img_array.shape[3]):\n"," masked_input[0,:,:,k] *= act_map_normalized\n"," masked_input_list.append(masked_input)\n"," masked_input_array = np.concatenate(masked_input_list, axis=0)\n"," # 4. feed masked inputs into CNN model and softmax\n"," #pred_from_masked_input_array = softmax(model.predict(masked_input_array))\n"," #print(masked_input_array.shape)\n"," pred_from_masked_input_array = softmax(model.predict(masked_input_array))\n"," #pred_from_masked_input_array = model.predict(masked_input_array)\n"," # 5. define weight as the score of target class\n"," weights = pred_from_masked_input_array[:,cls]\n"," # 6. get final class discriminative localization map as linear weighted combination of all activation maps\n","\n"," act_map_array=np.maximum(0,act_map_array)\n"," cam = np.dot(act_map_array[0,:,:,:], weights)\n"," cam = np.maximum(0, cam) # Passing through ReLU\n"," cam /= np.max(cam) # scale 0 to 1.0\n"," \n"," return cam\n","\n","def softmax(x):\n"," f = np.exp(x)/np.sum(np.exp(x), axis = 1, keepdims = True)\n"," return f\n","\n"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["def GuidedGradCam(model, img_array, layer_name, class_index, eps=1e-8):\n","\n"," gradModel = tf.keras.models.Model(\n","\t\t\tinputs=[model.inputs],\n","\t\t\toutputs=[model.get_layer(layer_name).output,\n","\t\t\t\tmodel.output])\n"," \n"," with tf.GradientTape() as tape:\n","\t\t\t# cast the image tensor to a float-32 data type, pass the\n","\t\t\t# image through the gradient model, and grab the loss\n","\t\t\t# associated with the specific class index\n"," inputs = tf.cast(img_array, tf.float32)\n"," (convOutputs, predictions) = gradModel(inputs)\n"," loss = predictions[:, class_index]\n","\t\t# use automatic differentiation to compute the gradients\n"," grads = tape.gradient(loss, convOutputs)\n"," \n"," # compute the guided gradients\n"," castConvOutputs = tf.cast(convOutputs > 0, \"float32\")\n"," castGrads = tf.cast(grads > 0, \"float32\")\n"," guidedGrads = castConvOutputs * castGrads * grads\n","\t\t# the convolution and guided gradients have a batch dimension\n","\t\t# (which we don't need) so let's grab the volume itself and\n","\t\t# discard the batch\n"," convOutputs = convOutputs[0]\n"," guidedGrads = guidedGrads[0]\n"," # compute the average of the gradient values, and using them\n","\t\t# as weights, compute the ponderation of the filters with\n","\t\t# respect to the weights\n"," weights = tf.reduce_mean(guidedGrads, axis=(0, 1))\n"," cam = tf.reduce_sum(tf.multiply(weights, convOutputs), axis=-1)\n"," \n"," # grab the spatial dimensions of the input image and resize\n","\t\t# the output class activation map to match the input image\n","\t\t# dimensions\n"," (w, h) = (img_array.shape[2], img_array.shape[1])\n"," heatmap = cv2.resize(cam.numpy(), (w, h))\n","\t\t# normalize the heatmap such that all values lie in the range\n","\t\t# [0, 1], scale the resulting values to the range [0, 255],\n","\t\t# and then convert to an unsigned 8-bit integer\n"," numer = heatmap - np.min(heatmap)\n"," denom = (heatmap.max() - heatmap.min()) + eps\n"," heatmap = numer / denom\n"," return heatmap\n","\n","\n","def sigmoid(x, a, b, c):\n"," return c / (1 + np.exp(-a * (x-b)))"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["from PIL import Image\n","# Dictionary where images and their probability are stored\n","images_dict={}\n","# vector to check for repeated images\n","repetead=[]\n","# Images that get a probability greater than 0.8 are in a file that follows the format:\n","# image name: probability: original image width: original image height\n","with open('name.txt') as f:\n"," while True:\n"," line = f.readline()\n"," line=line.rstrip()\n"," if not line:\n"," break\n"," image, auc, w, h=line.split(\":\")\n"," # The file is checked for repeated images (sanity comprobation)\n"," # and a probability is set (example of images above 0.8)\n"," probability=0.8\n"," if image not in repetead and float(auc)>probability:\n"," repetead.append(image)\n"," images_dict[image]=(float(auc), (w, h)) \n"," break \n"," else:\n"," # Show the repetead image\n"," print(image)"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["def get_IoU(truth_coords, pred_coords):\n"," # coords of intersection rectangle\n"," x1 = max(truth_coords[0], pred_coords[0])\n"," y1 = max(truth_coords[1], pred_coords[1])\n"," x2 = min(truth_coords[2], pred_coords[2])\n"," y2 = min(truth_coords[3], pred_coords[3])\n"," # area of intersection rectangle\n"," interArea = max(0, x2 - x1 + 1) * max(0, y2 - y1 + 1)\n"," # area of prediction and truth rectangles\n"," boxTruthArea = (truth_coords[2] - truth_coords[0] + 1) * (truth_coords[3] - truth_coords[1] + 1)\n"," boxPredArea = (pred_coords[2] - pred_coords[0] + 1) * (pred_coords[3] - pred_coords[1] + 1)\n"," # intersection over union \n"," iou=interArea /float(boxTruthArea + boxPredArea - interArea)\n"," return iou"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["# Function that normalizes de image\n","def normalize_image(array):\n","\tarray = tf.cast(array,tf.float32) # Convert to float\n","\tarray = tf.math.multiply(array, 1./255) # Convert to range 0-1\n","\tarray = tf.math.subtract(array, 0.449) # Normalize((x-mean)/std)\n","\tarray = tf.math.multiply(array, 1./0.226)\n","\treturn array\n","\n","# Trained model from which the heatmaps are to be obtained\n","model=tf.keras.models.load_model(\"model.h5\")\n","\n","heatmaps_guided_grad_cam=[]\n","heatmaps_grad_cam_plus=[]\n","heatmaps_score_cam=[]\n","heatmaps_faster_score_cam=[]\n","heatmaps_grad_cam=[]\n","heatmaps_dict={}\n","\n","for img_path in images_dict.keys():\n"," print(img_path)\n"," s,size=images_dict[img_path]\n"," _img=tf.io.decode_png(tf.io.read_file(img_path))\n"," # Resize image to (320,320)\n"," _img=tf.image.resize(_img, (320,320))\n"," # Normalize the image\n"," _img = normalize_image(_img)\n"," # Batch shape to model prediction\n"," _img = np.expand_dims(_img, axis=0)\n"," img = cv2.imread(img_path)\n","\n"," # Layer of the model from which the heatmaps are to be obtained \n"," layer=\"last layer model name\"\n"," \n"," \n"," # Output for which the heat map is to be obtained (disease)\n"," # For example if we train the a network which targets follow the next order: \n"," # Enlarged Cardiomediastinum,Cardiomegaly,Lung Opacity,Lung Lesion,Edema,Consolidation,\n"," #Pneumonia,Atelectasis,Pneumothorax,Pleural Effusion,Pleural Other,Fracture,Support Devices\n"," output=1\n","\n"," #Grad cam\n"," heatmap = xray_CAM(model,_img, layer, output)\n"," heatmap = cv2.resize(heatmap, (img.shape[1], img.shape[0]))\n"," heatmaps_grad_cam.append(heatmap)\n","\n"," plt.figure(figsize=(15, 15))\n"," ax = plt.subplot(1, 5, 1)\n"," plt.axis('off')\n"," impose=superimpose(img_path, heatmap)\n"," plt.imshow(impose)\n"," plt.show()\n"," \n","\n","\n"," #Guided Grad cam\n"," heatmap = GuidedGradCam(model,_img, layer, output)\n"," heatmap = cv2.resize(heatmap, (img.shape[1], img.shape[0]))\n"," heatmaps_guided_grad_cam.append(heatmap)\n","\n"," plt.figure(figsize=(15, 15))\n"," ax = plt.subplot(1, 5, 1)\n"," plt.axis('off')\n"," impose=superimpose(img_path, heatmap)\n"," plt.imshow(impose)\n"," plt.show()\n"," \n","\n"," #Grad cam plus\n"," heatmap = grad_cam_plus(model,_img, layer, output)\n"," heatmap = cv2.resize(heatmap, (img.shape[1], img.shape[0]))\n"," heatmaps_grad_cam_plus.append(heatmap)\n","\n"," \n"," plt.figure(figsize=(15, 15))\n"," ax = plt.subplot(1, 5, 1)\n"," plt.axis('off')\n"," impose=superimpose(img_path, heatmap)\n"," plt.imshow(impose)\n"," plt.show()\n","\n"," \n"," #Score cam\n"," heatmap = ScoreCam(model,_img, layer, output)\n"," heatmap = cv2.resize(heatmap, (img.shape[1], img.shape[0]))\n"," heatmaps_score_cam.append(heatmap)\n","\n","\n"," plt.figure(figsize=(15, 15))\n"," ax = plt.subplot(1, 5, 1)\n"," plt.axis('off')\n"," impose=superimpose(img_path, heatmap)\n"," plt.imshow(impose)\n"," plt.show()\n","\n"," #Faster score cam\n"," heatmap = ScoreCam(model,_img, layer, output, 10)\n"," heatmap = cv2.resize(heatmap, (img.shape[1], img.shape[0]))\n"," faster=heatmap\n"," heatmaps_faster_score_cam.append(heatmap)\n","\n"," \n"," plt.figure(figsize=(15, 15))\n"," ax = plt.subplot(1, 5, 1)\n"," plt.axis('off')\n"," impose=superimpose(img_path, heatmap)\n"," plt.imshow(impose)\n"," plt.show()\n","\n"," "]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["# Obtain Bounding boxes form heatmaps given a threshold\n","\n","def form_bboxes(heatmap, threshold):\n"," \n"," grey_heatmap = cv2.normalize(heatmap, None, alpha=0, beta=255, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_8U)\n"," ret,thresh = cv2.threshold(grey_heatmap,threshold,255,cv2.THRESH_BINARY)\n"," contours,hierarchy = cv2.findContours(thresh, cv2.RETR_LIST, cv2.CHAIN_APPROX_NONE)\n","\n"," # Value f the maximum pixel in the grey scale heatmap \n"," max=np.max(grey_heatmap)\n","\n"," not_box=True\n"," bbox_coords=[]\n"," global_maximum=0\n","\n"," \n"," for item in range(len(contours)):\n"," cnt = contours[item]\n"," x,y,w,h = cv2.boundingRect(cnt) # x, y is the top left corner, and w, h are the width and height respectively\n"," # coordinates in the image are upside down\n"," matriz=np.array(grey_heatmap[y:y+h, x:x+w])\n"," \n"," num_of_maximum_pixels=len(matriz[matriz==max])\n"," if num_of_maximum_pixels>=global_maximum:\n"," bbox_coords=[x,y,w,h]\n"," not_box=False\n"," global_maximum=num_of_maximum_pixels\n"," \n","\n"," # # If no contours are obtained with the given threshold, an empty list is returned. \n"," if not_box:\n"," return []\n"," else:\n"," return bbox_coords"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["new_df=pd.read_csv('bbox_nih.csv')\n","a=new_df[new_df[\"Finding Label\"]==\"Cardiomegaly\"]"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["# Construct dict with the bounding boxes form the heatmaps of grad-cam, guided grad-cam, grad-cam++, score-cam and faster score-cam\n","bbox_dict={}\n","for index, path in enumerate(images_dict.keys()): \n"," \n"," img = cv2.imread(path)\n"," bbox_dict[path]=[]\n","\n"," for heatmaps_list in [heatmaps_grad_cam, heatmaps_guided_grad_cam, heatmaps_grad_cam_plus, heatmaps_score_cam, heatmaps_faster_score_cam]:\n"," heatmap =heatmaps_list[index]\n"," # Heatmap cnormalization to apply threshold_multiotsu function\n"," h = cv2.normalize(heatmap, None, alpha=0, beta=255, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_8U)\n"," \n"," \n"," \n"," # Apply multiotusu algorith to obtain the theshols used to segment the heatmap to obtain the bounding box\n"," t0, t1, t2=threshold_multiotsu(h, classes=4)\n"," bbox = form_bboxes(heatmap, t2)\n"," \n"," [xp, yp, wp, hp]= bbox\n"," bbox_dict[path].append([xp, yp, wp, hp])\n"," "]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["# The NIH's csv where the annotations of the radiologists are\n","df=pd.read_csv('bbox_nih.csv')\n","#('bbox nih csv')\n","# Fix a pathology (example Cardiomegaly)\n","pathology_bounding_boxes=df[df[\"Finding Label\"]==\"Cardiomegaly\"]"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["import seaborn as snNew\n","import pandas as pdNew\n","import matplotlib.pyplot as pltNew\n","\n","total_iou=[]\n","\n","vector_heatmaps=[heatmaps_grad_cam, heatmaps_guided_grad_cam, heatmaps_grad_cam_plus, heatmaps_score_cam, heatmaps_faster_score_cam]\n","# len of matrix + 1 (Compare the heatmaps of the methods of the vector with the radiologist bbox too)\n","num=len(vector_heatmaps)+1\n","for key in (images_dict.keys()):\n"," matrix_iou=np.zeros((num,num))\n"," np.fill_diagonal(matrix_iou, 1)\n","\n"," img = cv2.imread(key)\n","\n"," csv_path=key.split(\"/\")[1]\n"," # radiologist coordinates of bounding box\n"," xr=int(np.round(pathology_bounding_boxes[pathology_bounding_boxes[\"Image Index\"]==csv_path][\"Bbox [x\"].tolist()[0],0))\n"," yr=int(np.round(pathology_bounding_boxes[pathology_bounding_boxes[\"Image Index\"]==csv_path][\"y\"].tolist()[0],0))\n"," wr=int(np.round(pathology_bounding_boxes[pathology_bounding_boxes[\"Image Index\"]==csv_path][\"w\"].tolist()[0],0))\n"," hr=int(np.round(pathology_bounding_boxes[pathology_bounding_boxes[\"Image Index\"]==csv_path][\"h]\"].tolist()[0],0))\n"," \n","\n"," for e, elem in enumerate(vector_heatmaps):\n"," \n"," [xb,yb,wb,hb]=bbox_dict[key][e]\n","\n"," iou=get_IoU([xr,yr,wr+xr,hr+yr], [xb, yb, xb+wb, yb+hb])\n"," e=e+1\n"," matrix_iou[e][0]=iou\n"," matrix_iou[0][e]=iou\n"," \n","\n"," for i in range(e,len(vector_heatmaps)):\n"," [xn,yn,wn,hn]=bbox_dict[key][i]\n"," i=i+1\n"," iou=get_IoU([xn,yn,wn+xn,hn+yn], [xb, yb, xb+wb, yb+hb])\n"," matrix_iou[e][i]=iou\n"," matrix_iou[i][e]=iou\n"," \n"," \n"," total_iou.append(matrix_iou)\n"," "]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["fig=plt.figure(figsize=(15,5))\n","snNew.set(font_scale=1.5)\n"," \n","ax = plt.subplot(1, 1, 1)\n","DetaFrame_cm = pdNew.DataFrame(np.round(np.sum(total_iou, axis=0)/len(total_iou),2), [\"R\", \"GC\", \"GGC\", \"G++\", \"SC\", \"FSC\"], [\"R\", \"GC\", \"GGC\", \"G++\", \"SC\", \"FSC\"])\n","snNew.heatmap(DetaFrame_cm, annot=True, ax=ax)\n","ax.set_title('IOU', fontsize=20)\n","\n","\n","plt.tight_layout()\n","plt.show()\n"]}],"metadata":{"colab":{"provenance":[]},"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.10.0"},"vscode":{"interpreter":{"hash":"916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1"}}},"nbformat":4,"nbformat_minor":0}
