{"cells":[{"metadata":{},"cell_type":"markdown","source":"For a discussion of the methods used here see [this](https://www.kaggle.com/gautham11/generating-scan-and-segemented-scan-jpg-fastai/notebook) kernel and if you want to understand what is going on you will have to look there. Note that the method does not work perfectly as can be seen by several dark images in the plots generated while making tfrecords."},{"metadata":{"trusted":true,"_kg_hide-output":true,"_kg_hide-input":true},"cell_type":"code","source":"!conda install -c conda-forge pillow -y\n!conda install -c conda-forge pydicom -y\n!conda install -c conda-forge gdcm -y\n!pip install pylibjpeg pylibjpeg-libjpeg","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"import os\nimport cv2\n\nimport pydicom\nimport pandas as pd\nimport numpy as np \nimport torch\nimport tensorflow as tf \nimport matplotlib.pyplot as plt \nfrom pathlib import Path\nimport scipy.ndimage as ndimage\nfrom skimage import measure, morphology, segmentation\nfrom scipy.ndimage.interpolation import zoom\nfrom PIL import Image \n\nfrom tqdm.notebook import tqdm\n%matplotlib inline","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"!pip install ../input/fastai2-wheels/fastcore-0.1.18-py3-none-any.whl -q\n!pip install ../input/fastai2-wheels/fastai2-0.0.17-py3-none-any.whl -q","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"from fastai2.basics import *\nfrom fastai2.medical.imaging import *\nimport cv2\nimport torch","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"data_dir = Path('/kaggle/input/osic-pulmonary-fibrosis-progression/')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"## Code from https://www.kaggle.com/aadhavvignesh/lung-segmentation-by-marker-controlled-watershed \n\ndef generate_markers(image):\n    \"\"\"\n    Generates markers for a given image.\n    \n    Parameters: image\n    \n    Returns: Internal Marker, External Marker, Watershed Marker\n    \"\"\"\n    \n    #Creation of the internal Marker\n    marker_internal = image < -400\n    marker_internal = segmentation.clear_border(marker_internal)\n    marker_internal_labels = measure.label(marker_internal)\n    \n    areas = [r.area for r in measure.regionprops(marker_internal_labels)]\n    areas.sort()\n    \n    if len(areas) > 2:\n        for region in measure.regionprops(marker_internal_labels):\n            if region.area < areas[-2]:\n                for coordinates in region.coords:                \n                       marker_internal_labels[coordinates[0], coordinates[1]] = 0\n    \n    marker_internal = marker_internal_labels > 0\n    \n    # Creation of the External Marker\n    external_a = ndimage.binary_dilation(marker_internal, iterations=10)\n    external_b = ndimage.binary_dilation(marker_internal, iterations=55)\n    marker_external = external_b ^ external_a\n    \n    # Creation of the Watershed Marker\n    marker_watershed = np.zeros(image.shape, dtype=np.int)\n    marker_watershed += marker_internal * 255\n    marker_watershed += marker_external * 128\n    \n    return marker_internal, marker_external, marker_watershed\n\n\ndef seperate_lungs(image, iterations = 1):\n    \"\"\"\n    Segments lungs using various techniques.\n    \n    Parameters: image (Scan image), iterations (more iterations, more accurate mask)\n    \n    Returns: \n        - Segmented Lung\n        - Lung Filter\n        - Outline Lung\n        - Watershed Lung\n        - Sobel Gradient\n    \"\"\"\n    \n    # Store the start time\n    # start = time.time()\n    \n    marker_internal, marker_external, marker_watershed = generate_markers(image)\n    \n    \n    '''\n    Creation of Sobel Gradient\n    '''\n    \n    # Sobel-Gradient\n    sobel_filtered_dx = ndimage.sobel(image, 1)\n    sobel_filtered_dy = ndimage.sobel(image, 0)\n    sobel_gradient = np.hypot(sobel_filtered_dx, sobel_filtered_dy)\n    sobel_gradient *= 255.0 / np.max(sobel_gradient)\n    \n    \n    '''\n    Using the watershed algorithm\n    \n    \n    We pass the image convoluted by sobel operator and the watershed marker\n    to morphology.watershed and get a matrix matrix labeled using the \n    watershed segmentation algorithm.\n    '''\n    watershed = morphology.watershed(sobel_gradient, marker_watershed)\n    \n    '''\n    Reducing the image to outlines after Watershed algorithm\n    '''\n    outline = ndimage.morphological_gradient(watershed, size=(3,3))\n    outline = outline.astype(bool)\n    \n    \n    '''\n    Black Top-hat Morphology:\n    \n    The black top hat of an image is defined as its morphological closing\n    minus the original image. This operation returns the dark spots of the\n    image that are smaller than the structuring element. Note that dark \n    spots in the original image are bright spots after the black top hat.\n    '''\n    \n    # Structuring element used for the filter\n    blackhat_struct = [[0, 0, 1, 1, 1, 0, 0],\n                       [0, 1, 1, 1, 1, 1, 0],\n                       [1, 1, 1, 1, 1, 1, 1],\n                       [1, 1, 1, 1, 1, 1, 1],\n                       [1, 1, 1, 1, 1, 1, 1],\n                       [0, 1, 1, 1, 1, 1, 0],\n                       [0, 0, 1, 1, 1, 0, 0]]\n    \n    blackhat_struct = ndimage.iterate_structure(blackhat_struct, iterations)\n    \n    # Perform Black Top-hat filter\n    outline += ndimage.black_tophat(outline, structure=blackhat_struct)\n    \n    '''\n    Generate lung filter using internal marker and outline.\n    '''\n    lungfilter = np.bitwise_or(marker_internal, outline)\n    lungfilter = ndimage.morphology.binary_closing(lungfilter, structure=np.ones((5,5)), iterations=3)\n    \n    '''\n    Segment lung using lungfilter and the image.\n    '''\n    segmented = np.where(lungfilter == 1, image, -2000*np.ones(image.shape))\n    \n    return segmented, lungfilter, outline, watershed, sobel_gradient\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"#DICOM Read utils\ndef fix_pxrepr(dcm):\n    if dcm.PixelRepresentation != 0 or dcm.RescaleIntercept<-100:\n        return dcm\n    x = dcm.pixel_array + 1000\n    px_mode = 4096\n    x[x>=px_mode] = x[x>=px_mode] - px_mode\n    dcm.PixelData = x.tobytes()\n    dcm.RescaleIntercept = -1000\n    return dcm\n\ndef read_dcm(path):\n    dcm = fix_pxrepr(Path(path).dcmread())\n    if dcm.Rows != 512 or dcm.Columns != 512: \n        dcm.zoom_to((512,512))\n    return dcm","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"bins = torch.tensor([-4096., -3024., -2048., -2000., -1109., -1025., -1024., -1023., -1020.,\n        -1017., -1014., -1011., -1008., -1005., -1003., -1001.,  -998.,  -996.,\n         -993.,  -991.,  -988.,  -985.,  -981.,  -978.,  -973.,  -969.,  -964.,\n         -958.,  -952.,  -946.,  -940.,  -933.,  -926.,  -919.,  -912.,  -904.,\n         -895.,  -886.,  -875.,  -862.,  -847.,  -829.,  -806.,  -777.,  -738.,\n         -684.,  -607.,  -501.,  -386.,  -301.,  -247.,  -210.,  -183.,  -163.,\n         -148.,  -136.,  -126.,  -117.,  -109.,  -101.,   -93.,   -85.,   -77.,\n          -69.,   -60.,   -51.,   -42.,   -33.,   -24.,   -15.,    -6.,     0.,\n            5.,    13.,    20.,    27.,    34.,    41.,    48.,    55.,    62.,\n           70.,    79.,    90.,   102.,   118.,   137.,   162.,   193.,   237.,\n          301.,   417.,   750.,  1278.])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"masked_bins = torch.tensor([-2048., -1000.,  -942.,  -908.,  -879.,  -847.,  -808.,  -751.,  -656.,\n         -467.,   -81.,   251.])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Generate tfrecords\nThis was helped by [this](https://www.kaggle.com/cdeotte/how-to-create-tfrecords) kernel."},{"metadata":{"trusted":true},"cell_type":"code","source":"groups = np.load('../input/cv-splits/groups.npy',allow_pickle=True)\n\n# Load meta data\nINPUT_FOLDER = '../input/osic-pulmonary-fibrosis-progression/train/'\ntrain = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/train.csv')\n# patients = os.listdir(INPUT_FOLDER)\npatients = train['Patient'].unique()\n\n# Make a dictionary that can read the number of images for a given patient\nimport os\ndirector = \"../input/osic-pulmonary-fibrosis-progression/train\"\nnum_im_dict = {}\nfor k, pid in enumerate(patients):\n    num_im_dict[pid] = len(os.listdir( director + \"/\" + pid ))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def _bytes_feature(value):\n  \"\"\"Returns a bytes_list from a string / byte.\"\"\"\n  if isinstance(value, type(tf.constant(0))):\n    value = value.numpy() # BytesList won't unpack a string from an EagerTensor.\n  return tf.train.Feature(bytes_list=tf.train.BytesList(value=[value]))\n\n# def _float_feature(value):\n#   \"\"\"Returns a float_list from a float / double.\"\"\"\n#   return tf.train.Feature(float_list=tf.train.FloatList(value=[value]))\n\ndef _floats_feature(value):\n  return tf.train.Feature(float_list=tf.train.FloatList(value=value.reshape(-1)))\n\ndef _int64_feature(value):\n  \"\"\"Returns an int64_list from a bool / enum / int / uint.\"\"\"\n  return tf.train.Feature(int64_list=tf.train.Int64List(value=[value]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def serialize_example(feature0, feature1, feature2, feature3):\n  feature = {\n      'image': _bytes_feature(feature0),\n      'image_name': _bytes_feature(feature1),\n      'image_dim': _int64_feature(feature2),\n      'slice': _int64_feature(feature3)\n  }\n  example_proto = tf.train.Example(features=tf.train.Features(feature=feature))\n  return example_proto.SerializeToString()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def process_image(path):\n    dcom = fix_pxrepr(Path(path).dcmread())\n    img = dcom.scaled_px.numpy()\n    if img.shape[0] != img.shape[1]:\n        mid_height = int(np.floor(img.shape[0]/2))\n        mid_width = int(np.floor(img.shape[1]/2))\n        img = img[(mid_height-256):(mid_height+256),(mid_width-256):(mid_width+256)]\n    if img.shape[0] != 512:\n        resize_factor = np.array([512,512]) / np.array(img.shape)\n        img = scipy.ndimage.interpolation.zoom(img, resize_factor, mode='nearest')\n\n    scaled_img = torch.tensor(img).hist_scaled(bins).numpy()\n\n    _, lungmask, _, _, _ = seperate_lungs(img)\n    masked_lung = torch.tensor(np.where(lungmask, img, -2048)).hist_scaled(masked_bins).numpy()\n    chest_lung = torch.tensor(np.where(~lungmask, img, -2048)).hist_scaled(masked_bins).numpy()\n\n    return np.stack((scaled_img,masked_lung,chest_lung))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(44,4,figsize=(20,44*5))\npax = ax.flatten()\n\nCT = len(groups)\nz=0; scaled_img = np.stack((np.eye(3),np.eye(3)))\nfor j,patients in enumerate(groups):\n    print(); print('Writing TFRecord %i of %i...'%(j,CT))\n    CT2 = len(patients)\n    tot_im_num = np.sum([num_im_dict[pid] for pid in patients])\n    if scaled_img.shape[1]!=scaled_img.shape[2]:print('Image size not equal');break\n    with tf.io.TFRecordWriter('train%.2i-%i-%i.tfrec'%(j,CT2,tot_im_num)) as writer:\n        for k,pt in enumerate(patients):\n            nim = num_im_dict[pt]\n            indx = int(nim/2)\n            name = str.encode(pt)\n            slc = 0; imcnt = 0\n            # If the image is not there don't count it, but do increase the slice number\n            while imcnt < nim-1:\n                try:\n                    slc+=1\n                    scaled_img = process_image(data_dir/'train/{}/{}.dcm'.format(pt,slc+1))\n                    imcnt+=1\n                    dimg = scaled_img.shape[1]\n                    img = cv2.imencode('.jpg', np.rollaxis(scaled_img,0,3)*255, (cv2.IMWRITE_JPEG_QUALITY, 94))[1].tostring()\n                    example = serialize_example(\n                        img, \n                        name,\n                        dimg,\n                        slc)\n                    writer.write(example)\n                    if slc == indx:\n                        pax[z].imshow(scaled_img[0], cmap=plt.cm.gray)\n                        z+=1\n                except FileNotFoundError: pass\n\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Verify TFRecords\nWe will verify the TFRecords we just made by using code from the Flower Comp starter notebook [here][1] to display the TFRecords below.\n\n[1]: https://www.kaggle.com/mgornergoogle/getting-started-with-100-flowers-on-tpu"},{"metadata":{"trusted":true},"cell_type":"code","source":"# numpy and matplotlib defaults\nnp.set_printoptions(threshold=15, linewidth=80)\nCLASSES = [0,1]\n\ndef batch_to_numpy_images_and_labels(data):\n    images, labels = data\n    numpy_images = images.numpy()\n    numpy_labels = labels.numpy()\n    #if numpy_labels.dtype == object: # binary string in this case, these are image ID strings\n    #    numpy_labels = [None for _ in enumerate(numpy_images)]\n    # If no labels, only image IDs, return None for labels (this is the case for test data)\n    return numpy_images, numpy_labels\n\ndef title_from_label_and_target(label, correct_label):\n    if correct_label is None:\n        return CLASSES[label], True\n    correct = (label == correct_label)\n    return \"{} [{}{}{}]\".format(CLASSES[label], 'OK' if correct else 'NO', u\"\\u2192\" if not correct else '',\n                                CLASSES[correct_label] if not correct else ''), correct\n\ndef display_one_flower(image, title, subplot, red=False, titlesize=16):\n    plt.subplot(*subplot)\n    plt.axis('off')\n    plt.imshow(image)\n    if len(title) > 0:\n        plt.title(title, fontsize=int(titlesize) if not red else int(titlesize/1.2), color='red' if red else 'black', fontdict={'verticalalignment':'center'}, pad=int(titlesize/1.5))\n    return (subplot[0], subplot[1], subplot[2]+1)\n    \ndef display_batch_of_images(databatch, predictions=None):\n    \"\"\"This will work with:\n    display_batch_of_images(images)\n    display_batch_of_images(images, predictions)\n    display_batch_of_images((images, labels))\n    display_batch_of_images((images, labels), predictions)\n    \"\"\"\n    # data\n    images, labels = batch_to_numpy_images_and_labels(databatch)\n    if labels is None:\n        labels = [None for _ in enumerate(images)]\n        \n    # auto-squaring: this will drop data that does not fit into square or square-ish rectangle\n    rows = int(math.sqrt(len(images)))\n    cols = len(images)//rows\n        \n    # size and spacing\n    FIGSIZE = 13.0\n    SPACING = 0.1\n    subplot=(rows,cols,1)\n    if rows < cols:\n        plt.figure(figsize=(FIGSIZE,FIGSIZE/cols*rows))\n    else:\n        plt.figure(figsize=(FIGSIZE/rows*cols,FIGSIZE))\n    \n    # display\n    for i, (image, label) in enumerate(zip(images[:rows*cols], labels[:rows*cols])):\n        title = label\n        correct = True\n        if predictions is not None:\n            title, correct = title_from_label_and_target(predictions[i], label)\n        dynamic_titlesize = FIGSIZE*SPACING/max(rows,cols)*40+3 # magic formula tested to work from 1x1 to 10x10 images\n        subplot = display_one_flower(image, title, subplot, not correct, titlesize=dynamic_titlesize)\n    \n    #layout\n    plt.tight_layout()\n    if label is None and predictions is None:\n        plt.subplots_adjust(wspace=0, hspace=0)\n    else:\n        plt.subplots_adjust(wspace=SPACING, hspace=SPACING)\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def decode_image(image_data):\n    image = tf.image.decode_jpeg(image_data, channels=3)\n    image = tf.cast(image, tf.float32) / 255.0  # convert image to floats in [0, 1] range\n    image = tf.reshape(image, [*IMAGE_SIZE, 3]) # explicit size needed for TPU\n    return image\n\ndef read_labeled_tfrecord(example):\n    LABELED_TFREC_FORMAT = {\n        \"image\": tf.io.FixedLenFeature([], tf.string), # tf.string means bytestring\n        \"image_name\": tf.io.FixedLenFeature([], tf.string),  # shape [] means single element\n    }\n    example = tf.io.parse_single_example(example, LABELED_TFREC_FORMAT)\n    image = decode_image(example['image'])\n    label = example['image_name']\n    return image, label # returns a dataset of (image, label) pairs\n\ndef load_dataset(filenames, labeled=True, ordered=False):\n    # Read from TFRecords. For optimal performance, reading from multiple files at once and\n    # disregarding data order. Order does not matter since we will be shuffling the data anyway.\n\n    ignore_order = tf.data.Options()\n    if not ordered:\n        ignore_order.experimental_deterministic = False # disable order, increase speed\n\n    dataset = tf.data.TFRecordDataset(filenames, num_parallel_reads=AUTO) # automatically interleaves reads from multiple files\n    dataset = dataset.with_options(ignore_order) # uses data as soon as it streams in, rather than in its original order\n    dataset = dataset.map(read_labeled_tfrecord)\n    # returns a dataset of (image, label) pairs if labeled=True or (image, id) pairs if labeled=False\n    return dataset\n\ndef get_training_dataset():\n    dataset = load_dataset(TRAINING_FILENAMES, labeled=True)\n    dataset = dataset.repeat() # the training dataset must repeat for several epochs\n    dataset = dataset.shuffle(2048)\n    dataset = dataset.batch(BATCH_SIZE)\n    dataset = dataset.prefetch(AUTO) # prefetch next batch while training (autotune prefetch buffer size)\n    return dataset\n\ndef count_data_items(filenames):\n    # the number of data items is written in the name of the .tfrec files, i.e. flowers00-230.tfrec = 230 data items\n    n = [int(re.compile(r\"-([0-9]*)\\.\").search(filename).group(1)) for filename in filenames]\n    return np.sum(n)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# INITIALIZE VARIABLES\nIMAGE_SIZE= [512,512]; BATCH_SIZE = 32\nAUTO = tf.data.experimental.AUTOTUNE\nTRAINING_FILENAMES = tf.io.gfile.glob('train*.tfrec')\nprint('There are %i train images'%count_data_items(TRAINING_FILENAMES))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# DISPLAY TRAIN IMAGES\ntraining_dataset = get_training_dataset()\ntraining_dataset = training_dataset.unbatch().batch(20)\ntrain_batch = iter(training_dataset)\n\ndisplay_batch_of_images(next(train_batch))","execution_count":null,"outputs":[]}],"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}