extracting_patches_from_cbis_ddsm_1.py

# -*- coding: utf-8 -*-
"""Extracting_Patches_from_CBIS_DDSM_1.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1s-2TRdN0BBYqbU4WEbLLfvasgM96BWb_
"""

from google.colab import drive
drive.mount('/content/drive')

!pip install SimpleITK==2.1.1
!pip install pydicom

import pydicom
import SimpleITK as sitk
import numpy as np
import pandas as pd
import cv2
from tqdm.notebook import tqdm
import copy
import os
import glob
import shutil
import random
from sklearn.utils import shuffle
import matplotlib.pyplot as plt

def select_largest_obj(img_bin, lab_val=255, fill_holes=False, 
                       smooth_boundary=False, kernel_size=15):
    '''Select the largest object from a binary image and optionally
    fill holes inside it and smooth its boundary.
    Args:
        img_bin(2D array): 2D numpy array of binary image.
        lab_val([int]): integer value used for the label of the largest 
                        object. Default is 255.
        fill_holes([boolean]): whether fill the holes inside the largest 
                               object or not. Default is false.
        smooth_boundary([boolean]): whether smooth the boundary of the 
                                    largest object using morphological 
                                    opening or not. Default is false.
        kernel_size([int]): the size of the kernel used for morphological 
                            operation.
    '''
    n_labels, img_labeled, lab_stats, _ = cv2.connectedComponentsWithStats(
        img_bin, connectivity=8, ltype=cv2.CV_32S)
    largest_obj_lab = np.argmax(lab_stats[1:, 4]) + 1
    largest_mask = np.zeros(img_bin.shape, dtype=np.uint8)
    largest_mask[img_labeled == largest_obj_lab] = lab_val
    if fill_holes:
        bkg_locs = np.where(img_labeled == 0)
        bkg_seed = (bkg_locs[0][0], bkg_locs[1][0])
        img_floodfill = largest_mask.copy()
        h_, w_ = largest_mask.shape
        mask_ = np.zeros((h_ + 2, w_ + 2), dtype=np.uint8)
        cv2.floodFill(img_floodfill, mask_, seedPoint=bkg_seed, newVal=lab_val)
        holes_mask = cv2.bitwise_not(img_floodfill)  # mask of the holes.
        largest_mask = largest_mask + holes_mask
    if smooth_boundary:
        kernel_ = np.ones((kernel_size, kernel_size), dtype=np.uint8)
        largest_mask = cv2.morphologyEx(largest_mask, cv2.MORPH_OPEN, kernel_)
    
    cnts, _ = cv2.findContours(largest_mask.astype(np.uint8), cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
    cnt = max(cnts, key = cv2.contourArea)
    x, y, w, h = cv2.boundingRect(cnt)
    return x, y, w, h, largest_mask

def cropp_image(img, roi, threshold, truncate = False):
    _, binary_img = cv2.threshold(img, threshold, 
                                    maxval=255, type=cv2.THRESH_BINARY)
    x, y, w, h, mask_img = select_largest_obj(binary_img, lab_val=255, 
                                       fill_holes=False, 
                                       smooth_boundary=True, kernel_size=2)
    croped_img, br_mask = img[y:y+h, x:x+w], mask_img[y:y+h, x:x+w]
    croped_roi = roi[y:y+h, x:x+w]

    if truncate:
        Pmin = np.percentile(croped_img[br_mask!=0], 5)
        Pmax = np.percentile(croped_img[br_mask!=0], 99)
        truncated = np.clip(croped_img,Pmin, Pmax)  
        normalized = (truncated - Pmin)/(Pmax - Pmin)
        normalized[br_mask==0]=0
        filtered_img = normalized
    else:
        filtered_img = croped_img
    filtered_img = ((filtered_img - filtered_img.min()) / (filtered_img.max() - filtered_img.min())*255).astype(np.uint8)
    croped_roi = (croped_roi>0).astype(np.uint8)
    return filtered_img, croped_roi

import numpy as np
from scipy.ndimage import zoom


def clipped_zoom(img, mask, zoom_factor, **kwargs):

    h, w = img.shape[:2]

    zoom_tuple = (zoom_factor,) * 2 + (1,) * (img.ndim - 2)
    zh = int(np.round(h / zoom_factor))
    zw = int(np.round(w))
    top = (h - zh) // 2
    left = (w - zw) // 2

    out_img = zoom(img[top:top+zh, left:left+zw], zoom_tuple, **kwargs)
    out_mask = zoom(mask[top:top+zh, left:left+zw], zoom_tuple, **kwargs)
    return out_img, out_mask

def read_img_mask(img_path, mask_path):
    img = sitk.ReadImage(img_path)
    img_hu = sitk.GetArrayFromImage(img)[0]
    img = (((img_hu - img_hu.min())/(img_hu.max() - img_hu.min()))*255).astype(np.uint8)

    mask = sitk.ReadImage(mask_path) 
    mask = sitk.GetArrayFromImage(mask)[0]
    mask = (mask > 0).astype(np.uint8)

    img, mask = clipped_zoom(img, mask, 1.2)
    img, mask = cropp_image(img, mask, 5)
    return img, mask

root_path = '/content/drive/MyDrive/Breast Cancer Datasets/'

"""# Create ROI MetaData"""

mass1 = pd.read_csv(root_path+'CBIS/CBIS-DDSM-All-doiJNLP-zzWs5zfZ/mass_case_description_train_set.csv')
mass2 = pd.read_csv(root_path+'CBIS/CBIS-DDSM-All-doiJNLP-zzWs5zfZ/mass_case_description_test_set.csv')
mass = mass1.append(mass2, ignore_index=True)
mass.head()

calc1 = pd.read_csv(root_path+'CBIS/CBIS-DDSM-All-doiJNLP-zzWs5zfZ/calc_case_description_train_set.csv')
calc2 = pd.read_csv(root_path+'CBIS/CBIS-DDSM-All-doiJNLP-zzWs5zfZ/calc_case_description_test_set.csv')
calc = calc1.append(calc2, ignore_index=True)
calc.head()

cbis = mass.append(calc, ignore_index=True)
cbis.head()

meta = pd.read_csv(root_path+'CBIS/CBIS-DDSM-All-doiJNLP-zzWs5zfZ/metadata.csv')
meta_full = meta.loc[meta['Series Description']=='full mammogram images']
meta_roi = meta.loc[meta['Series Description']=='ROI mask images']
print(len(meta_full))
print(len(meta_roi))

meta_roi.head()

patients = cbis['patient_id'].values
print(len(patients))
patients = np.unique(patients)
print(len(patients))

patient_id = []
side = []
cc_path = []
mlo_path = []
cc_mask = []
mlo_mask = []
pathology = []
density = []
assessment = []
ab_type = []
for i in tqdm(range(len(patients))):
    df_patient = cbis.loc[cbis['patient_id']==patients[i]]
    for s in ['LEFT', 'RIGHT']:
        df_side = df_patient.loc[df_patient['left or right breast'] == s]
        if len(np.unique(df_side['image view'].values)) >= 2:
            try:
                cc = df_side.loc[df_side['image view']=='CC']['image file path'].values[0]
                cc = glob.glob(root_path+'CBIS/CBIS-DDSM-All-doiJNLP-zzWs5zfZ/'+
                                meta_full.loc[meta_full['Subject ID']==cc.split('/')[0]]['File Location'].values[0][2:]+'/*.dcm')[0]
                cc = cc.replace(root_path, '')

                cc_roi = df_side.loc[df_side['image view']=='CC']['ROI mask file path'].values[0]
                cc_roi = glob.glob(root_path+'CBIS/CBIS-DDSM-All-doiJNLP-zzWs5zfZ/'+
                                    meta_roi.loc[meta_roi['Subject ID']==cc_roi.split('/')[0]]['File Location'].values[0][2:]+'/*.dcm')
                cc_roi_des = []
                for ds in cc_roi:
                    dss =pydicom.dcmread(ds, stop_before_pixels=True)
                    try:
                        if dss.SeriesDescription == 'ROI mask images':
                            cc_roi_des.append(True)
                        else:
                            cc_roi_des.append(False)
                    except:
                        cc_roi_des.append(True)
                cc_roi = [cc_roi[nnn] for nnn in range(len(cc_roi)) if cc_roi_des[nnn]]
                cc_roi = [xx.replace(root_path, '') for xx in cc_roi]
                cc_roi = cc_roi[0]

                mlo = df_side.loc[df_side['image view']=='MLO']['image file path'].values[0]
                mlo = glob.glob(root_path+'CBIS/CBIS-DDSM-All-doiJNLP-zzWs5zfZ/'+
                                meta_full.loc[meta_full['Subject ID']==mlo.split('/')[0]]['File Location'].values[0][2:]+'/*.dcm')[0]
                mlo = mlo.replace(root_path, '')

                mlo_roi = df_side.loc[df_side['image view']=='MLO']['ROI mask file path'].values[0]
                mlo_roi = glob.glob(root_path+'CBIS/CBIS-DDSM-All-doiJNLP-zzWs5zfZ/'+
                                    meta_roi.loc[meta_roi['Subject ID']==mlo_roi.split('/')[0]]['File Location'].values[0][2:]+'/*.dcm')
                mlo_roi_des = []
                for ds in mlo_roi:
                    dss =pydicom.dcmread(ds, stop_before_pixels=True)
                    try:
                        if dss.SeriesDescription == 'ROI mask images':
                            mlo_roi_des.append(True)
                        else:
                            mlo_roi_des.append(False)
                    except:
                        mlo_roi_des.append(True)
                mlo_roi = [mlo_roi[nnn] for nnn in range(len(mlo_roi)) if mlo_roi_des[nnn]]
                mlo_roi = [xx.replace(root_path, '') for xx in mlo_roi]
                mlo_roi = mlo_roi[0]
            except:
                continue
            cc_path.append(cc)
            cc_mask.append(cc_roi)
            mlo_path.append(mlo)
            mlo_mask.append(mlo_roi)
            patient_id.append(patients[i])
            side.append(s)
            density.append(df_side['breast_density'].values[0])
            pathh = df_side['pathology'].values[0]
            if '_' in pathh:
                pathh = 'BENIGN'
            pathology.append(pathh.lower().capitalize())
            assessment.append(df_side['assessment'].values[0])
            ab_type.append(df_side['abnormality type'].values[0])

index = np.arange(len(patient_id))
columns = ['PatientID', 'Side', 'CC', 'CC-Mask', 'MLO', 'MLO-Mask', 'Pathology', 'Density', 'BI-RADS', 'Type']
df = pd.DataFrame(index=index, columns=columns)

df['PatientID'] = patient_id
df['Side'] = side
df['CC'] = cc_path
df['CC-Mask'] = cc_mask
df['MLO'] = mlo_path
df['MLO-Mask'] = mlo_mask
df['Pathology'] = pathology
df['Density'] = density
df['BI-RADS'] = assessment
df['Type'] = ab_type

df

clinical = df

n = 1000
sample = clinical.iloc[n]
CC = sample['CC']
CC_mask = sample['CC-Mask']
MLO = sample['MLO']
MLO_mask = sample['MLO-Mask']
abnormality = sample['Type']
label = sample['Pathology']

print(sample['PatientID'])
print(label, abnormality)
print('Images: ')

dicoms = [CC, MLO]
for i in range(len(dicoms)):
    img = sitk.ReadImage(root_path+'/'+dicoms[i])
    img = sitk.GetArrayFromImage(img)[0]
    if i == 0:
        masks = CC_mask
    else:
        masks = MLO_mask

    plt.figure(figsize=(8, 4))
    plt.subplot(1, 2, 1)
    plt.imshow(img, cmap='gray')
    plt.axis('off')
    plt.title('original')
    roi = sitk.ReadImage(root_path+'/'+masks)
    roi = sitk.GetArrayFromImage(roi)[0]
    roi = (roi > 0).astype(np.uint8)
    plt.subplot(1, 2, 2)
    plt.imshow(img, cmap='gray')
    plt.imshow(roi, alpha=0.5, cmap='gray')
    plt.axis('off')
    plt.show()

clinical.to_csv(root_path+'CBIS/CBIS-DDSM-All-doiJNLP-zzWs5zfZ/cbis-roi-metadata.csv')

"""# Read ROIs"""

clinical = pd.read_csv(root_path+'CBIS/CBIS-DDSM-All-doiJNLP-zzWs5zfZ/cbis-roi-metadata.csv')
clinical.head()

img_path = list(clinical['CC'].values)
mask_path = list(clinical['CC-Mask'].values)

img_path.extend(list(clinical['MLO'].values))
mask_path.extend(list(clinical['MLO-Mask'].values))

print(len(img_path))
print(len(mask_path))

n = 100
img, mask = read_img_mask(root_path+img_path[n], root_path+mask_path[n])
plt.figure()
plt.subplot(1,2,1)
plt.imshow(img, cmap='gray')
plt.axis('off')
plt.subplot(1,2,2)
plt.imshow(img, cmap='gray')
plt.imshow(mask, alpha=0.5, cmap='gray')
plt.axis('off')
plt.show()

"""# Create Sample Patches"""

clinical = pd.read_csv(root_path+'CBIS/CBIS-DDSM-All-doiJNLP-zzWs5zfZ/cbis-roi-metadata.csv')
clinical.head()

def add_img_margins(img, margin_size):
    '''Add all zero margins to an image
    '''
    enlarged_img = np.zeros((int(img.shape[0]+margin_size*2), 
                             int(img.shape[1]+margin_size*2)))
    enlarged_img[int(margin_size):int(margin_size+img.shape[0]), 
                 int(margin_size):int(margin_size+img.shape[1])] = img
    return enlarged_img

def crop_val(v, minv, maxv):
    v = v if v >= minv else minv
    v = v if v <= maxv else maxv
    return v

def overlap_patch_roi(patch_center, patch_size, roi_mask, 
                      add_val=1000, cutoff=.5):
    x1,y1 = (patch_center[0] - patch_size/2, 
             patch_center[1] - patch_size/2)
    x2,y2 = (patch_center[0] + patch_size/2, 
             patch_center[1] + patch_size/2)
    x1 = crop_val(x1, 0, roi_mask.shape[1])
    y1 = crop_val(y1, 0, roi_mask.shape[0])
    x2 = crop_val(x2, 0, roi_mask.shape[1])
    y2 = crop_val(y2, 0, roi_mask.shape[0])
    roi_area = (roi_mask>0).sum()
    roi_patch_added = roi_mask.copy()
    roi_patch_added[int(y1):int(y2), int(x1):int(x2)] += add_val
    patch_area = (roi_patch_added>=add_val).sum()
    inter_area = (roi_patch_added>add_val).sum().astype(np.float32)
    return (inter_area/roi_area > cutoff or inter_area/patch_area > cutoff)

def sample_patches(img, roi_mask, patch_directory, patch_size=256,
                   pos_cutoff=.75, neg_cutoff=.35,
                   nb_bkg=10, nb_abn=10, start_sample_nb=0, verbose=False):

    img = add_img_margins(img, patch_size/2)
    roi_mask = add_img_margins(roi_mask, patch_size/2)
    # Get ROI bounding box.
    roi_mask_8u = roi_mask.astype(np.uint8)
    contours,_ = cv2.findContours(
        roi_mask_8u.copy(), cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
    cont_areas = [ cv2.contourArea(cont) for cont in contours ]
    idx = np.argmax(cont_areas)  # find the largest contour.
    rx,ry,rw,rh = cv2.boundingRect(contours[idx])
    if verbose:
        M = cv2.moments(contours[idx])
        try:
            cx = int(M['m10']/M['m00'])
            cy = int(M['m01']/M['m00'])
        except ZeroDivisionError:
            cx = rx + int(rw/2)
            cy = ry + int(rh/2)
    rng = np.random.RandomState(12345)
    # Sample abnormality first.
    sampled_abn = 0
    nb_try = 0
    os.mkdir(patch_directory+'/abnormalities/')
    while sampled_abn < nb_abn:
        if nb_abn > 1:
            x = rng.randint(rx, rx + rw)
            y = rng.randint(ry, ry + rh)
            nb_try += 1
            if nb_try >= 1000:
                print("Nb of trials reached maximum, decrease overlap cutoff by 0.05")
                pos_cutoff -= .05
                nb_try = 0
                if pos_cutoff <= .0:
                    raise Exception("overlap cutoff becomes non-positive, "
                                    "check roi mask input.")
        else:
            x = cx
            y = cy
        # import pdb; pdb.set_trace()
        if nb_abn == 1 or overlap_patch_roi((x,y), patch_size, roi_mask, 
                                            cutoff=pos_cutoff):
            patch = img[int(y - patch_size/2):int(y + patch_size/2), 
                        int(x - patch_size/2):int(x + patch_size/2)]
            patch = patch.astype(np.uint8)
            patch_name = patch_directory+'/abnormalities/'+'abn_patch_{}.jpg'.format(sampled_abn)
            cv2.imwrite(patch_name, patch)
            sampled_abn += 1
            nb_try = 0
    # Sample background.
    sampled_bkg = start_sample_nb
    os.mkdir(patch_directory+'/backgrounds/')
    while sampled_bkg < start_sample_nb + nb_bkg:
        x = rng.randint(patch_size/2, img.shape[1] - patch_size/2)
        y = rng.randint(patch_size/2, img.shape[0] - patch_size/2)
        if not overlap_patch_roi((x,y), patch_size, roi_mask, cutoff=neg_cutoff):
            patch = img[int(y - patch_size/2):int(y + patch_size/2), 
                        int(x - patch_size/2):int(x + patch_size/2)]
            patch = patch.astype(np.uint8)
            patch_name = patch_directory+'/backgrounds/'+'bkg_patch_{}.jpg'.format(sampled_bkg)
            cv2.imwrite(patch_name, patch)
            sampled_bkg += 1

p = clinical['PatientID'].values
p, n = np.unique(p, return_counts=True)
np.unique(n, return_counts=True)

os.mkdir('CBIS-Patches_10')

for i in tqdm(range(len(clinical))):
    p_id = clinical['PatientID'].values[i]
    side = clinical['Side'].values[i][0]

    p_folder = 'CBIS-Patches_10/{}/'.format(p_id)
    p_folders = os.listdir('CBIS-Patches_10/')
    if p_id not in p_folders:
        os.mkdir(p_folder)

    s_folder = p_folder+'{}/'.format(side)
    os.mkdir(s_folder)
    
    for view in ['CC', 'MLO']:
        view_folder = s_folder+'{}/'.format(view)
        os.mkdir(view_folder)
        img_path = root_path+clinical[view].values[i]
        mask_path = root_path+clinical['{}-Mask'.format(view)].values[i]
        patch_direc = view_folder
        img, mask = read_img_mask(img_path, mask_path)
        sample_patches(img, mask, patch_direc)

shutil.make_archive('/content/drive/MyDrive/Breast Cancer Preprocessed Datasets/CBIS-Patches_10',
                    'zip',
                    '/content/CBIS-Patches_10'
                    )

# os.mkdir('CBIS-Patches')

# !unzip -d CBIS-Patches/ '/content/drive/MyDrive/Breast Cancer Preprocessed Datasets/CBIS-Patches_1'
# !unzip -d CBIS-Patches/ '/content/drive/MyDrive/Breast Cancer Preprocessed Datasets/CBIS-Patches_2'
# !unzip -d CBIS-Patches/ '/content/drive/MyDrive/Breast Cancer Preprocessed Datasets/CBIS-Patches_3'
# !unzip -d CBIS-Patches/ '/content/drive/MyDrive/Breast Cancer Preprocessed Datasets/CBIS-Patches_4'
# !unzip -d CBIS-Patches/ '/content/drive/MyDrive/Breast Cancer Preprocessed Datasets/CBIS-Patches_5'
# !unzip -d CBIS-Patches/ '/content/drive/MyDrive/Breast Cancer Preprocessed Datasets/CBIS-Patches_6'

# patients = glob.glob('CBIS-Patches/*')
# print(len(patients))

# patientID = []
# Side = []
# PatchFolder_abn = []
# PatchFolder_bck = []
# Label = []


# for i in tqdm(range(len(patients))):
#     sides = os.listdir(patients[i])

#     for s in sides:
#         if s == 'L':
#             s_t = 'LEFT'
#         else:
#             s_t = 'RIGHT'
#         dff = clinical.loc[(clinical['PatientID']==patients[i].split('/')[-1]) & (clinical['Side']==s_t)]
#         label = dff['Pathology'].values[0] + ' ' + dff['Type'].values[0]

#         patch_folders_abn = glob.glob(patients[i]+'/'+s+'/*/abnormalities/')
#         patch_folders_bck = glob.glob(patients[i]+'/'+s+'/*/backgrounds/')

#         patientID.append(patients[i].split('/')[-1])
#         Side.append(s)
#         PatchFolder_abn.append(patch_folders_abn)
#         PatchFolder_bck.append(patch_folders_bck)
#         Label.append(label)

# index = np.arange(len(patientID))
# columns = ['PatientID', 'Side', 'PatchFolder_abn', 'PatchFolder_bck', 'Label']

# df = pd.DataFrame(index=index, columns=columns)

# df['PatientID'] = patientID
# df['Side'] = Side
# df['PatchFolder_abn'] = PatchFolder_abn
# df['PatchFolder_bck'] = PatchFolder_bck
# df['Label'] = Label
# df

# from sklearn.model_selection import train_test_split

# x = list(df.index)
# y = df['Label'].values

# x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2, random_state=2)

# Training_Tag = np.empty(len(x), dtype='object')

# for i in x_train:
#     Training_Tag[i] = 'Train'

# for i in x_test:
#     Training_Tag[i] = 'Evaluation'

# Training_Tag

# df['Training_Tag'] = Training_Tag
# df

# df.to_csv('/content/drive/MyDrive/Breast Cancer Preprocessed Datasets/cbis-patch-labels.csv', index=False)

# shutil.make_archive('/content/drive/MyDrive/Breast Cancer Preprocessed Datasets/CBIS-Patches',
#                     'zip',
#                     '/content/CBIS-Patches'
#                     )