{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Part 1 - Tradition Image processing for classification: feature extraction\n\nThe notebook is an educational experiment rather than a competitive one. The notebook explores traditional image processing methods for extracting features from images to use in classification. I've only used a subset of available techniques here, several other features can also be derived from the images.\n\nNeedless to mention, modern deep learning approaches will perform significantly better. It was fun to play around.","execution_count":null},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":true},"cell_type":"code","source":"import pandas as pd \nimport numpy as np\nimport matplotlib\nimport matplotlib.pyplot as plt\nfrom tqdm import tqdm_notebook\nfrom matplotlib.patches import Rectangle\nimport seaborn as sns\nimport pydicom as dcm\n%matplotlib inline \nplt.set_cmap(plt.cm.bone)\nIS_LOCAL = True\nimport os\nimport cv2\nimport skimage\nfrom skimage import feature, filters\nfrom tqdm import tqdm\n\nPATH=\"../input/rsna-pneumonia-detection-challenge\"\n\nprint(os.listdir(PATH))","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"class_info_df = pd.read_csv(PATH+'/stage_2_detailed_class_info.csv')\ntrain_labels_df = pd.read_csv(PATH+'/stage_2_train_labels.csv')\ntrain_class_df = train_labels_df.merge(class_info_df, left_on='patientId', right_on='patientId', how='inner')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_class_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Utility functions","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_images(data):\n imgs = []\n for path in data['patientId']:\n patientImage = path + '.dcm'\n imagePath = os.path.join(PATH,\"stage_2_train_images/\", patientImage)\n img = dcm.read_file(imagePath).pixel_array\n imgs.append(img)\n return imgs\n\ndef imshow_gray(img):\n plt.figure(figsize=(12,7))\n return plt.imshow(img, cmap='gray')\n \ndef imshow_with_labels(img, patient_id):\n rows = train_labels_df[train_labels_df['patientId'] == patient_id]\n for row in rows.itertuples(): \n x, y, w, h = row.x, row.y, row.width, row.height\n x, y, w, h = map(int, [x,y,w,h])\n cv2.rectangle(img, (x,y), (x+w,y+h), 255, 2)\n plt.figure(figsize=(12,7))\n return plt.imshow(img, cmap='gray')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Sample a few images and process them","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"test_df = train_class_df[train_class_df['Target']==1].sample(4)\nbox = test_df.loc[test_df.index, ['x', 'y', 'width', 'height']]\ntest = load_images(test_df[0:3])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"idx = 1\nimg = test[idx]\nimshow_with_labels(img.copy(), test_df.iloc[idx,0])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Image enhancement","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### Histogram equalization","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"equ = cv2.equalizeHist(test[idx])\nax = imshow_gray(equ)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Equalization presents a good contrast of the lungs, and further accents the presence of opacity","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### Image sharpening","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# apply hpf kernel\nhpf_kernel = np.full((3, 3), -1)\nhpf_kernel[1,1] = 9\nim_hp = cv2.filter2D(equ, -1, hpf_kernel)\n\n# use unsharpen mask filter\nim_us = skimage.filters.unsharp_mask(equ)\n\nfig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16,9))\nax1.imshow(im_hp, cmap='gray')\nax1.set_title('High pass filter')\nax2.imshow(im_us, cmap='gray')\nax2.set_title('unsharpen mask filter')\nfig.suptitle('Image sharpening')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"High pass show more detail, while unmask is still a little blurry, not ideal for feature extraction","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### Thresholding\n","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# otsu thresholding\nret, otsu = cv2.threshold(cv2.GaussianBlur(im_hp,(7,7),0),0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)\n# local thresholding\nlocal = im_hp > skimage.filters.threshold_local(im_hp, 5)\nmean = im_hp > skimage.filters.threshold_mean(im_hp)\n\nfig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(16,9))\nax1.imshow(otsu, cmap='gray')\nax1.set_title('Otsu thresholding')\nax2.imshow(local, cmap='gray')\nax2.set_title('Local thresholding')\nfig.suptitle('Image sharpening')\nax3.imshow(mean, cmap='gray')\nax3.set_title('Mean thresholding')\nfig.suptitle('Image thresholding')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Otsu thresholding helps extract smooth edges, so we choose it and it segments the lungs portion better","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### Edge detection","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"sobel = filters.sobel(otsu)\ncanny = feature.canny(otsu/255)\n\nfig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16,9))\nax1.imshow(canny, cmap='gray')\nax1.set_title('Canny edge detection')\nax2.imshow(sobel, cmap='gray')\nax2.set_title('Sobel operator')\nfig.suptitle('Edge detection')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Sobel fiter extracts the edges better than canny","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### Lung segmentation\n\nFind and pull out contours","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"contours, hier = cv2.findContours((sobel * 255).astype('uint8'),cv2.RETR_EXTERNAL,cv2.CHAIN_APPROX_SIMPLE)\nprint('Contours found ', len(contours))\n\nsrt_contours = sorted(contours, key=lambda x: x.shape[0], reverse=True)\nselect_contour = srt_contours[0] # probably not the best assumption\n\ntest = img.copy()\nimg_contour = cv2.drawContours(test, [select_contour], 0, (255,0,0), thickness=3)\n\nimshow_gray(img_contour)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"After identifying the lung segment we can extract the center of moment of this segment. Since, all the X-ray images are from the same dimension, this can be a valid feature for prediction","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"M = cv2.moments(select_contour)\ncx = int(M['m10'] / M['m00'])\ncy = int(M[\"m01\"] / M[\"m00\"])\n\ntest = img.copy()\ncv2.circle(test, (cx, cy), 7, (255, 255, 255), -1)\nimshow_gray(test)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"But there are images were the subject is slightly rotated, and differently sized so we want center of moment to be invariant to rotation and scale. So we pick **Hu moments**. We will also log the moments to make it easy to compare and drop the 3rd moment as it depends on the other moments and 7th moment as it distinguishes mirror images and there are no flipped images in the dataset","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_hu_moments(contour):\n M = cv2.moments(select_contour)\n hu = cv2.HuMoments(M).ravel().tolist()\n del hu[2]\n del hu[-1]\n log_hu = [-np.sign(a)*np.log10(np.abs(a)) for a in hu]\n return log_hu\n\nget_hu_moments(select_contour)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Feature extraction\nWe derive the following features from our above results to build a classifier for pneumonia detection:\n* Area of opacity\n* Perimeter of visible lung regions\n* Irregularity index\n* Equivalent diameter\n* Mean, sd of unenhanced image\n* Hu moments\n","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def area(img):\n # binarized image as input\n return np.count_nonzero(img)\n\ndef perimeter(img):\n # edges of the image as input\n return np.count_nonzero(img)\n\ndef irregularity(area, perimeter):\n # area and perimeter of the image as input, also called compactness\n I = (4 * np.pi * area) / (perimeter ** 2)\n return I\n\ndef equiv_diam(area):\n # area of image as input\n ed = np.sqrt((4 * area) / np.pi)\n return ed\n\ndef get_hu_moments(contour):\n # hu moments except 3rd and 7th (5 values)\n M = cv2.moments(contour)\n hu = cv2.HuMoments(M).ravel().tolist()\n del hu[2]\n del hu[-1]\n log_hu = [-np.sign(a)*np.log10(np.abs(a)) for a in hu]\n return log_hu","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### The image preprocessing and feature extraction pipeline","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def extract_features(img):\n mean = img.mean()\n std_dev = img.std()\n \n # hist equalization\n equalized = cv2.equalizeHist(img)\n \n # sharpening\n hpf_kernel = np.full((3, 3), -1)\n hpf_kernel[1,1] = 9\n sharpened = cv2.filter2D(equalized, -1, hpf_kernel)\n \n # thresholding\n ret, binarized = cv2.threshold(cv2.GaussianBlur(sharpened,(7,7),0),0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)\n \n # edge detection\n edges = skimage.filters.sobel(binarized)\n \n # moments from contours\n contours, hier = cv2.findContours((edges * 255).astype('uint8'),cv2.RETR_EXTERNAL,cv2.CHAIN_APPROX_SIMPLE)\n select_contour = sorted(contours, key=lambda x: x.shape[0], reverse=True)[0]\n \n # feature extraction\n ar = area(binarized)\n per = perimeter(edges)\n irreg = irregularity(ar, per)\n eq_diam = equiv_diam(ar)\n hu = get_hu_moments(select_contour)\n \n return (mean, std_dev, ar, per, irreg, eq_diam, *hu)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# test the function\nextract_features(img)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Save generated features","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Load the required data\nWe select only the normal and pneumonia images for model building","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"pneumonia_ids = train_labels_df[train_labels_df['Target'] == 1]['patientId'].unique()\npneumonia_labels = [1] * len(pneumonia_ids)\n\nnormal_ids = class_info_df[class_info_df['class'] == 'Normal']['patientId'].unique()\nnormal_labels = [0] * len(normal_ids)\n\ndata = dict()\ndata['patientId'] = np.concatenate((pneumonia_ids, normal_ids))\ndata['target'] = np.concatenate((pneumonia_labels, normal_labels))\n\nprint(f'Pneumonia images: {len(pneumonia_ids)}\\nNormal images: {len(normal_ids)}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Generate features\nFor each ID generate the features from the image and store it in the dataset","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"from tqdm import tqdm","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"features = []\n\nfor path in tqdm(data['patientId']):\n patientImage = path + '.dcm'\n imagePath = os.path.join(PATH,\"stage_2_train_images/\", patientImage)\n img = dcm.read_file(imagePath).pixel_array\n feats = extract_features(img)\n features.append(feats)\n\ndata['features'] = features","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df = pd.DataFrame(data)\ndf.to_csv('img_features.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"After the features are generated, they can be loaded and trained using machine learning models to perform classification. Results with different ML algorithms are presented in the the [2nd part of this notebook](https://www.kaggle.com/suryathiru/2-tradition-image-processing-model-training)","execution_count":null}],"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":4,"nbformat_minor":4}