<jupyter_start><jupyter_text>Downloading the SIIM-ACR Pneumothorax Segmentation datasetThe data is comprised of images in DICOM format and annotations in the form of image IDs and run-length-encoded (RLE) masks. Some of the images contain instances of pneumothorax (collapsed lung), which are indicated by encoded binary masks in the annotations. Some training images have multiple annotations.<jupyter_code>!kaggle datasets download -d jesperdramsch/siim-acr-pneumothorax-segmentation-data
!7z x lung-segmentation-for-siimacr-pneumothorax.zip
!7z x siim-acr-pneumothorax-segmentation-data.zip
import numpy as np
import pandas as pd
import tensorflow as tf
import cv2
from skimage import transform
from skimage.filters import threshold_otsu
from skimage.util import img_as_ubyte
from skimage import morphology, io, color, exposure, img_as_float, transform
from skimage.feature import local_binary_pattern
# To calculate a normalized histogram
from skimage.util.shape import view_as_blocks
from skimage.transform import rescale, resize
from scipy.stats import itemfreq
from sklearn.preprocessing import normalize
import pydicom
from glob import glob
from tqdm import tqdm<jupyter_output><empty_output><jupyter_text>Pre-processing(Lung Boundary Segmentation)  Lung Boundary Detection using UNet modelUNet is a Fully Convolutional Neural Network used for semantic segmentation.The model has 5 million parameters and has been trained on the Montgomery dataset.The model takes a 256x256 floating point image as input and returns a masking image of the same datatype and dimensions.Following cell depicts the structure in detail.<jupyter_code>model_dir = '/content/trained_model.hdf5'
lung_seg_model = tf.keras.models.load_model(model_dir, custom_objects=None, compile=True)
lung_seg_model.summary()<jupyter_output>WARNING:tensorflow:From /usr/local/lib/python3.6/dist-packages/tensorflow/python/ops/nn_impl.py:180: add_dispatch_support.<locals>.wrapper (from tensorflow.python.ops.array_ops) is deprecated and will be removed in a future version.
Instructions for updating:
Use tf.where in 2.0, which has the same broadcast rule as np.where
WARNING:tensorflow:Error in loading the saved optimizer state. As a result, your model is starting with a freshly initialized optimizer.
Model: "model_1"
__________________________________________________________________________________________________
Layer (type)                    Output Shape         Param #     Connected to                     
==================================================================================================
input_1 (InputLayer)            [(None, 256, 256, 1) 0                                            
__________________________________________________________________________________________________
conv2d_1 (Conv2D)       [...]<jupyter_text>The following function does these steps:1.   Resizes the image to (256,256)2.   Applies Histogram Equalization3.   Normalizes the image by subtracting the mean and dividing it by the std dev.<jupyter_code>def get_lung_seg_tensor(img, batch_size, seg_size, n_channels):
    
        'Generates data containing batch_size samples' # X : (n_samples, *dim, n_channels)
        # Initialization
        X = np.empty((batch_size, seg_size, seg_size, n_channels))

        # Process Image
        #pixel_array = pydicom.read_file(file_path).pixel_array
        pixel_array = img
        image_resized = cv2.resize(pixel_array, (seg_size, seg_size))
        image_resized = exposure.equalize_hist(image_resized)
        image_resized = np.array(image_resized, dtype=np.float64)
        image_resized -= image_resized.mean()
        image_resized /= image_resized.std()
        # Store Image
        X[0,] = np.expand_dims(image_resized, axis=2)

        return X<jupyter_output><empty_output><jupyter_text>The following function takes an image and a mask image, returns with only the Region of Interest(ROI) present based on the OTSU threshold of the mask image.In our case the ROI are the lung lobes in the Chest X-ray.<jupyter_code>def crop_lungs(org_img,thresh_img):
  im_size = org_img.shape
  if(org_img.shape != thresh_img):
    thresh_img = transform.resize(thresh_img,im_size)
  otsu_thresh = threshold_otsu(thresh_img)
  new_im = np.full(im_size,255)
  #print(org_img.shape," ",thresh_img.shape)
  for i in range(im_size[0]):
    for j in range(im_size[1]):
      if(thresh_img[i][j] > otsu_thresh):
        new_im[i][j] = org_img[i][j]
  #new_im = org_img>thresh_img
  return new_im
seg_size = 256
test_img = "/content/dicom-images-train/1.2.276.0.7230010.3.1.2.8323329.5682.1517875189.580785/1.2.276.0.7230010.3.1.3.8323329.5682.1517875189.580784/1.2.276.0.7230010.3.1.4.8323329.5682.1517875189.580786.dcm"
img = get_lung_seg_tensor(pydicom.read_file(test_img).pixel_array,1,seg_size,1)
seg_mask = lung_seg_model.predict(img).reshape((seg_size,seg_size))
from matplotlib import pyplot as plt
%matplotlib inline
fig, (ax1, ax2, ax3) = plt.subplots(1, 3,figsize = (15,15))
fig.suptitle('Figure 1.Localizing the ROI using the UNet model')
ax1.imshow(pydicom.read_file(test_img).pixel_array,cmap = 'bone')
ax1.title.set_text('(a)The original image')
ax2.imshow(seg_mask, cmap = 'bone')
ax2.title.set_text('(b)The Detected ROI')
ax3.title.set_text('(c)Extracted ROI using OTSU thresholding')
im3  = crop_lungs(pydicom.read_file(test_img).pixel_array,seg_mask)
ax3.imshow(im3, cmap = 'bone')
plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
from skimage import feature
class LBP:
  def __init__(self,points = 24, radius = 8):
    self.points = points
    self.radius = radius
  def feature_vector(self,image,eps = 1e-7,method = 'uniform'):
    lbp = feature.local_binary_pattern(image, self.points, self.radius, method= "uniform" )
    (hist, _) = np.histogram(lbp.ravel(),bins=np.arange(0, self.points + 3),range=(0, self.points + 2))
    hist = hist.astype("float")
    hist /= (hist.sum() + eps)
    return hist<jupyter_output><empty_output><jupyter_text>An example of the given annotated labelsThe annotated labels are given in the following formatIf 'EncodedPixels' = -1 then no Pneumothorax*   If 'EncodedPixels' = -1 then no Pneumothorax*   Else 'EncodedPixels' = RLE(Pneumothorax)In the following example ImageId 0 and 1 do not have pneumothorax whereas ImageIdf 2,3 and 3 have pneumothorax<jupyter_code>df = pd.read_csv('train-rle.csv')
df.head()
print("Total no. of examples in the SIIM-ACR Pneumothorax Segmentation dataset are ",len(df.index))
import sys
sys.path.insert(0, '/content/dicom-images-train')
import glob
import os
all_paths = []
all_paths = sorted(glob.glob('/content/dicom-images-train/**/*.dcm',recursive = True))<jupyter_output><empty_output><jupyter_text>using 8 neighboring points with 1 pixel radius for unifrom LBP<jupyter_code>seg_size = 256
eps = 1e-7
ULBP = LBP(points = 8 , radius = 1)
training_counter = 0<jupyter_output><empty_output><jupyter_text>Feature ExtractionIn machine learning, Feature Extraction is the process of arriving at a set of derived values that are intended to be informative and non-redundant from an inital set of raw values.The end product of this process is called a *feature vector*.The *feature vector* needs to be such that the classifier gets all the information it needs in a manner in which it understands. Local Binary PatternsLocal Binary Pattern (LBP) is a simple yet very efficient texture operator which labels the pixels of an image by thresholding the neighborhood of each pixel and considers the result as a binary number.LBP can be used to extract texture features from an image in the following way1.   Divide the examined image into cells2.   For each pixel in a cell, compare the pixel to each of its 8 neighbors 3.   Follow the pixels along a circle, i.e. clockwise4.   Where the center pixel's value is greater than the neighbor's value, write "0". Otherwise, write "1". This gives an 8-digit binary number5.   Compute the histogram, over the cell, of the frequency of each "number" occurring6. Normalize the histogram7. Concatenate (normalized) histograms of all cells. This gives a feature vector for the entire image.![alt text](https://encrypted-tbn0.gstatic.com/images?q=tbn%3AANd9GcTo-iaSOU_-hWBBLAWFeu6KIFPUcEuyXks2FM_4Pe05ETPaMLP3)**Figure 1**.Example of LBP computation for the pixel with value *'4'* results in the binary string *'001110100'* Uniform Local Binary PatternsA Local Binary Pattern is said to be [*uniform*](https://www.hindawi.com/journals/isrn/2013/429347/) if and only if has at most two 1->0 or 0->1 transitions when traversed **circularly** Concept of uniform patterns1.   Greatly reduces the size of the feature vector when it comes to higher number of sampling points2.   [Experiments](http://www.scholarpedia.org/article/Local_Binary_Patterns) have shown that uniform patterns occur more frequently in texture images than non-uniform patterns<jupyter_code>train_X = []
train_Y = []
block_size = 64
block_step = block_size
img_size = 512
training_counter = 0
for path in all_paths:
  #if (training_counter%500 == 0):
  #  print("So far length of Xtrain ",len(train_X)," length of Ytrain ",len(train_Y))
  img = pydicom.read_file(path).pixel_array
  img_norm = get_lung_seg_tensor(img,1,seg_size,1)
  seg_mask = lung_seg_model.predict(img_norm).reshape((seg_size,seg_size))
  seg_mask = resize(seg_mask,output_shape=(img_size,img_size))
  cropped = crop_lungs(img,seg_mask)
  tiles = view_as_blocks(cropped, block_shape=(block_size, block_size))
  #tiles = get_tile_images(cropped, 128, 128)
  X_i = np.array([ULBP.feature_vector(j) for i in tiles for j in i]).flatten()
  train_X.append(X_i)
  
  y =  df.loc[df['ImageId'] == path.split('/')[-1][:-4]].values.tolist()
  try:
    if '-1' in y[0][1]:
      train_Y.append(0)
    else:
      train_Y.append(1)
  except IndexError:
    train_Y.append(0)
  training_counter = training_counter + 1
X = np.array(train_X)
Y = np.array(train_Y)
X.shape
Y.shape<jupyter_output><empty_output><jupyter_text>Splitting the Dataset1.   **Training set**—A subset to train a model(80% of the total dataset).2.   **Test set**—A subset to test the trained model(20% of the total dataset)![alt text](https://developers.google.com/machine-learning/crash-course/images/PartitionTwoSets.svg)**Figure 2**. Visualization of slicing a single data set into a training set and test set.<jupyter_code>from sklearn.model_selection import train_test_split
X_train,X_test,Y_train,Y_test= train_test_split(X,Y,test_size = 0.2)<jupyter_output><empty_output><jupyter_text>Classification using Support Vector MachineA Support Vector Machine is defined as a constrained optimization problem to maximize the margin of separating hyperplanes.>$ y_{i} =(x_{i}^{T}\beta+\beta_{0}) $The values of $\beta$ are optimized by the following constrained optimiazation problem>$\min_{\beta_{0},\beta_{1}}{||\beta||}$subject to>$y_{i}(x_{i}^{T}\beta+\beta_{0})\geq1- \epsilon_{i}$>$ \forall_{i} \in (1,N) $$\epsilon_{i}$ = the amount by which $(x_{i}^{T}\beta+\beta_{0})$ is on the wrong sideUsing implementation from Scikit-learn: Machine Learning in Python, Pedregosa et al., JMLR 12, pp. 2825-2830, 2011.  Training the SVM model<jupyter_code>from sklearn.svm import LinearSVC
model = LinearSVC(C=100.0, random_state=42)
model.fit(X_train,Y_train)<jupyter_output>/usr/local/lib/python3.6/dist-packages/sklearn/svm/base.py:929: ConvergenceWarning: Liblinear failed to converge, increase the number of iterations.
  "the number of iterations.", ConvergenceWarning)<jupyter_text>Evaluating the modelAccuracy and Senesitivity(True Positive Rate) can be calculated from the following equations:>$Accuracy = \frac{TP+TN}{TP+TN+FP+FN}$>$Sensitivity =\frac{TP}{TP+FN} $<jupyter_code>Y_pred = model.predict(X_test)
from sklearn.metrics import accuracy_score
from sklearn.metrics import confusion_matrix
tn, fp, fn, tp = confusion_matrix(Y_test, Y_pred).ravel()
sensitivity = tp/(tp+fn)
print("Accuracy is (ULBP Feature Extraction) with Lung segmentation with SVM classifier",accuracy_score(Y_test,Y_pred))
result_dict = {'Y_true':Y_test,'Y_pred':Y_pred}
result_df = pd.DataFrame(result_dict)
result_df.to_csv('Pneumo_ULBP_SVM_Lung_Seg.csv')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10607426838265431, 'fraction_numerical': 0.028690876301336175, 'mean_word_length': 5.227667177131189, 'pattern_counts': {'":': 0, '<': 40, '<?xml version=': 0, '>': 50, 'https://': 3, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 3, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '209719', 'license': 'MIT', 'n_tokens_mistral': 4177, 'n_tokens_neox': 3890, 'n_words': 1129, 'path': 'Pneumothorax_Classification_with_ULBP,SVM&Lung_Segmentation.ipynb', 'repo_name': 'aaryapatel007/Pneumothorax-Classification-and-Segmentation', 'size': '12199'}