{"cells":[{"metadata":{},"cell_type":"markdown","source":"## Patient Class\nAn object oriented way of loading a DICOM image, each patient containing the following\n* A 4D image: $\\mathbb{R}^d$ where the dimensionality $d$ is given as $slices \\times time \\times height \\times width$\n* Pixel to millimeter ratio (width- and height-wise)\n* Spacing between slices (depth-wise)\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"import pydicom as dicom\nimport re\nimport os\nimport numpy as np\n\nclass Patient(object):\n    def __init__(self, directory, subdir):\n        # deal with any intervening directories\n        while True:\n            subdirs = next(os.walk(directory))[1]\n            if len(subdirs) == 1:\n                directory = os.path.join(directory, subdirs[0])\n            else:\n                break\n\n        slices = []\n        for s in subdirs:\n            m = re.match(\"sax_(\\d+)\", s)\n            if m is not None:\n                slices.append(int(m.group(1)))\n\n        slices_map = {}\n        first = True\n        times = []\n        for s in slices:\n            files = next(os.walk(os.path.join(directory, \"sax_%d\" % s)))[2]\n            offset = None\n\n            for f in files:\n                m = re.match(\"IM-(\\d{4,})-(\\d{4})\\.dcm\", f)\n                if m is not None:\n                    if first:\n                        times.append(int(m.group(2)))\n                    if offset is None:\n                        offset = int(m.group(1))\n\n            first = False\n            slices_map[s] = offset\n\n        self.directory = directory\n        self.time = sorted(times)\n        self.slices = sorted(slices)\n        self.slices_map = slices_map\n        self.name = subdir\n\n    def _filename(self, s, t):\n        fname = os.path.join(self.directory,\n                                 \"sax_%d\" % s, \n                                 \"IM-%04d-%04d.dcm\" % (self.slices_map[s], t))\n        return fname\n\n    def _read_dicom_image(self, filename):\n        d = dicom.read_file(filename)\n        img = d.pixel_array\n        return np.array(img)\n\n    def _read_all_dicom_images(self):\n        f1 = self._filename(self.slices[0], self.time[0])\n        f2 = self._filename(self.slices[1], self.time[0])\n        \n        d1 = dicom.read_file(f1)\n        d2 = dicom.read_file(f2)\n        \n        (x, y) = d1.PixelSpacing\n        (x, y) = (float(x), float(y))\n        self.col_scaling = x\n        self.row_scaling = y\n        \n        # try a couple of things to measure distance between slices\n        try:\n            dist = np.abs(d2.SliceLocation - d1.SliceLocation)\n        except AttributeError:\n            try:\n                dist = d1.SliceThickness\n            except AttributeError:\n                dist = 8  # better than nothing...\n\n        # 4D image array\n        self.images = np.array([[self._read_dicom_image(self._filename(d, i))\n                                for i in self.time]\n                                for d in self.slices])\n        \n        # Distance between slices in mm\n        self.dist = dist\n        \n        # Calculate depth as distance between slices times no. of slices\n        self.deph_mm = self.dist * (self.images.shape[0] - 1)\n        \n        # Area scaling, mm per pixel\n        self.area_multiplier = x * y\n        \n        # Orientation\n        self.orientation = d1.ImageOrientationPatient\n        \n    def load(self):\n        self._read_all_dicom_images()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Loading Patients\n* Load all patients using the Patient class above"},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_patient(patient_id, root_dir=None):\n    if not root_dir: \n        root_dir =  os.path.join('..', 'input', 'train', 'train')\n    patient_id = str(patient_id)\n    base_path = os.path.join(root_dir, patient_id)\n    try:\n        tdata = Patient(base_path, patient_id)\n        tdata.load()\n        # If data does not contain 4 dimensions, throw it away\n        if len(tdata.images.shape) == 4:\n            return tdata\n    except (ValueError, TypeError, IndexError, AttributeError, FileNotFoundError):\n        print('Patient %s could not be loaded.' % patient_id)\n        return None\n    \ndef load_multiple_patients(patient_ids=False, root_dir=None, verbose=False):\n    \"\"\"\n    :param patient_ids: ids of patients to load [list of integers]\n    :param root_dir: name of root dir, defaults to Kaggle root directory [string]\n    :param verbose: Whether to print every patient id when loading [boolean]\n    :return: list of [Patient] objects\n    \"\"\"\n    # If no ids are specified load all from 1-500\n    if not patient_ids:\n        patient_ids = range(1, 501)\n    patient_list = []\n    for pid in patient_ids:\n        if verbose:\n            print('Loading patient %i...' % pid)\n        p_data = load_patient(pid, root_dir=root_dir)\n        if p_data:\n            patient_list.append(p_data)\n    return patient_list","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Fourier Transform\n* Used for identifying regions with movement over time, in this case the left heart ventricle"},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\n\n# Based on https://gist.github.com/ajsander/fb2350535c737443c4e0#file-tutorial-md\ndef fourier_time_transform_slice(image_3d):\n    '''\n    3D array -> 2D array\n    [slice, height, width] -> [height, width]\n    Returns (width, height) matrix\n    Fourier transform for 3d data (time,height,weight)\n    '''\n    # Apply FFT to the selected slice\n    fft_img_2d = np.fft.fftn(image_3d)[1, :, :]\n    return np.abs(np.fft.ifftn(fft_img_2d))\n\n\ndef fourier_time_transform(patient_images):\n    '''\n    4D array -> 3D array (compresses time dimension)\n    Concretely, [slice, time, height, width] -> [slice, height, width]\n    Description: Fourier transform for analyzing movement over time.\n    '''\n\n    ftt_image = np.array([\n        fourier_time_transform_slice(patient_slice)\n        for patient_slice in patient_images\n    ])\n    return ftt_image","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Segmentation\n* Here, K means is used for finding a fitting threshhold for segmentation\n* The image is then segmented using the threshhold, effectively making every pixel either foreground (white = 1) or background (black = 0)\n* Lastly, by using an average of segmented pixel intesities, we identify the region of interest"},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\nfrom skimage.morphology import binary_dilation, binary_erosion, binary_opening, binary_closing, disk\nfrom skimage.filters import threshold_otsu\n\n\ndef thresh_segmentation(patient_img):\n    \"\"\"Returns matrix\n    Segmententation of patient_img with k-means\n    \"\"\"\n    #Z = np.float32(np.ravel(patient_img))\n    #criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 10, 1.0)\n    #flags = cv2.KMEANS_RANDOM_CENTERS\n    #compactness, labels, centers = cv2.kmeans(Z, 2, None, criteria, 10, flags)\n    #center = np.uint8(centers)\n    thresh = threshold_otsu(patient_img)\n    binary = patient_img > thresh\n    return binary\n\ndef segment_multiple(patient_img):\n    \"\"\"Returns list\n    List of segmented slices with function thresh_segmentation()\n    \"\"\"\n    num_slices, height, width = patient_img.shape\n    segmented_slices = np.zeros((num_slices, height, width))\n\n    for i in range(num_slices):\n        seg_slice = thresh_segmentation(patient_img[i])\n        if seg_slice.sum() > seg_slice.size * 0.5:\n            seg_slice = 1 - seg_slice\n        segmented_slices[i] = seg_slice\n\n    return segmented_slices\n\ndef roi_mean_yx(patient_img):\n    \"\"\"Returns mean(y) and mean(x) [double]\n    Mean coordinates in segmented patients slices.\n    This function performs erosion to get a better result.\n    Original: See https://nbviewer.jupyter.org/github/kmader/Quantitative-Big-Imaging-2019/blob/master/Lectures/06-ShapeAnalysis.ipynb\n    \"\"\"\n    seg_slices = segment_multiple(patient_img)\n    num_slices = seg_slices.shape[0]\n    y_all, x_all = np.zeros(num_slices), np.zeros(num_slices)\n    neighborhood = disk(2)\n    \n    for i,seg_slice in enumerate(seg_slices):\n        # Perform erosion to get rid of wrongly segmented small parts\n        seg_slices_eroded = binary_erosion(seg_slice, neighborhood) \n        \n        # Filter out background of slice, after erosion [background=0, foreground=1]\n        y_coord, x_coord = seg_slices_eroded.nonzero()\n        \n        # Save mean coordinates of foreground \n        y_all[i], x_all[i] = np.mean(y_coord), np.mean(x_coord)\n    \n    # Return mean of mean foregrounds - this gives an estimate of ROI coords.\n    mean_y = int(np.mean(y_all))\n    mean_x = int(np.mean(x_all))\n    return mean_y, mean_x","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Histogram Normalize\nApply histogram normalization to each 2d image in the 4d image\n\n* source: https://scikit-image.org/docs/dev/auto_examples/color_exposure/plot_equalize.html"},{"metadata":{"trusted":true},"cell_type":"code","source":"from skimage import exposure\n\ndef histogram_normalize_4d(images, clip_limit=0.03):\n    slices, time, _, _ = images.shape\n    norm_imgs_4d = np.empty(images.shape)\n    for i in range(slices):\n        for j in range(time):\n            norm_imgs_4d[i,j] = exposure.equalize_adapthist(images[i,j].astype(np.uint16), \n                                                            clip_limit=clip_limit)\n    return norm_imgs_4d","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Rescale Patient Images\nPatient data has been gathered on different devices, resulting in different image dimensions across patients. However, all DICOM images contain metadata about the scaling of the images, which we will use to normalize patient images.\nNext, we would like to remove unnecessary data, i.e. everything that is not the heart, since this cuts down on the input size for the data analysis.\n\nThe pre-processing is therefore a 2-step process:\n* Rescale patient images, such that 1 pixel = 1 mm\n* Crop out Region of Interest (Heart)"},{"metadata":{"trusted":true},"cell_type":"code","source":"import cv2\n\ndef rescale_patient_4d_imgs(patient):\n    img_4d = patient.images\n    if len(img_4d.shape) < 4: raise Exception(\"Patient images are not 4D!\")\n    num_slices, time, _, _ = img_4d.shape\n    \n    # Extract scaled DICOM width/height multipliers\n    # http://dicom.nema.org/dicom/2013/output/chtml/part03/sect_10.7.html\n    fx, fy = patient.col_scaling, patient.row_scaling\n    \n    # Rescale the first 2d image, in order to find out the resulting dimensions\n    example_img = cv2.resize(src=img_4d[0,0], dsize=None, fx=fx, fy=fy)\n    scaled_height, scaled_width = example_img.shape\n    scaled_imgs = np.zeros((num_slices, time, scaled_height, scaled_width))\n    \n    for i in range(num_slices):\n        for j in range(time):\n            scaled_imgs[i,j] = cv2.resize(src=img_4d[i,j], dsize=None, fx=fx, fy=fy)\n    \n    return scaled_imgs\n\ndef crop_roi(img, dim_y, dim_x, cy, cx):\n    \"\"\"\n    Crops an image from the given coords (cy, cx), such that the resulting img is of\n    dimensions [dim_y, dim_x], i.e. height and width.\n    Resulting image is filled out from top-left corner, and remaining pixels are left black.\n    \"\"\"\n    cy, cx = int(round(cy)), int(round(cx))\n    h, w = img.shape\n    if dim_x > w or dim_y > h: raise ValueError('Crop dimensions larger than image dimension!')\n    new_img = np.zeros((dim_y, dim_x))\n    dx, dy = int(dim_x / 2), int(dim_y / 2)\n    dx_odd, dy_odd = int(dim_x % 2 == 1), int(dim_y % 2 == 1)\n\n    # Find boundaries for cropping [original img]\n    dx_left = max(0, cx - dx)\n    dx_right = min(w, cx + dx + dx_odd)\n    dy_up = max(0, cy - dy)\n    dy_down = min(h, cy + dy + dy_odd)\n\n    # Find how many pixels to fill out in new image\n    range_x = dx_right - dx_left\n    range_y = dy_down - dy_up\n    \n\n    # Fill out new image from top left corner\n    # Leave pixels outside range as 0's (black)\n    new_img[0:range_y, 0:range_x] = img[dy_up:dy_down, dx_left:dx_right]\n    return new_img\n\ndef crop_heart(images_4d, heart_pixel_size=200):\n    # Find center for cropping\n    ft_imges = fourier_time_transform(images_4d)\n    y, x = roi_mean_yx(ft_imges)\n    \n    # Create new 4d image array\n    num_slices, time, h, w = images_4d.shape\n    heart_cropped_img_4d = np.zeros((num_slices, time, heart_pixel_size, heart_pixel_size))\n    \n    for i in range(num_slices):\n        for j in range(time):\n            heart_cropped_img_4d[i,j] = crop_roi(images_4d[i,j], heart_pixel_size, heart_pixel_size, y, x)\n    \n    return heart_cropped_img_4d\n\ndef rotate_images_210_deg(images_4d, orientation):\n    \"\"\"\n    Return 4d image\n    Params 4d numpy, int\n    Idea from: kaggle.com/c/second-annual-data-science-bowl/discussion/19378\n    Description: \n                Rotates image if orientation angle is -30 degreees, which ensures\n                that the left ventricle is in the top right corner of the image.\n    \"\"\"\n    angle = np.arctan2(orientation[:3], orientation[:3]) / np.pi * 180 - 75\n    rotation_needed = angle[2] > (-210)\n    \n    # Check if rotation needed\n    if rotation_needed:\n        # Calculate resulting dimensions for numpy array\n        slices, time, _, _ = images_4d.shape\n        rot_width, rot_height = np.rot90(images_4d[0,0], k=1).shape\n        rot_images = np.zeros((slices, time, rot_width, rot_height))\n        \n        # Rotate images\n        for i in range(slices):\n            for j in range(time):\n                rot_images[i,j] = np.rot90(images_4d[i,j], k=1)\n        return rot_images\n    \n    # Otherwise if no rotation needed, return original images\n    return images_4d","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"## Segmenting the Left Ventricle"},{"metadata":{"trusted":true},"cell_type":"code","source":"from skimage.morphology import opening, disk\nfrom scipy.ndimage import distance_transform_edt\nfrom skimage.morphology import watershed\nfrom skimage.feature import peak_local_max\n\n# Code from: https://nbviewer.jupyter.org/github/kmader/Quantitative-Big-Imaging-2019/blob/master/Lectures/07-ComplexObjects.ipynb\ndef watershed_img(image):\n    # Distance map\n    image_dmap = distance_transform_edt(image)\n    # Distance peaks\n    image_peaks = label(peak_local_max(image_dmap, indices=False, footprint=np.ones((40, 40)),labels=image, exclude_border=True))\n    # Watershed first once\n    ws_labels = watershed(-image_dmap, image_peaks, mask=image)\n    \n    # Reomve small segments\n    label_area_dict = {i: np.sum(ws_labels == i)for i in np.unique(ws_labels[ws_labels > 0])}\n    clean_label_maxi = image_peaks.copy()\n    lab_areas = list(label_area_dict.values())\n    # Remove 20 percentile\n    area_cutoff = np.percentile(lab_areas, 15)\n    for i, k in label_area_dict.items():\n        if k <= area_cutoff:\n            clean_label_maxi[clean_label_maxi == i] = 0\n    # Watershed again\n    ws_labels = watershed(-image_dmap, clean_label_maxi, mask=image)\n\n    return ws_labels","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from skimage.measure import label\n\ndef labeled_segmented_images(images):\n    \"\"\"\n    Returns numpy array (4d)\n    Segments image and used watershed for labeling.\n    \"\"\"\n    \n    num_slices, time, height, width = images.shape\n    segmented_slices = np.zeros((num_slices, time, height, width))\n    \n    # Iterate over all slices and whole timeseries for images\n    for i in range(num_slices):\n        for j in range(time):\n            # Segmentation\n            seg_slice = thresh_segmentation(images[i,j])\n            \n            # Makes all segmented images same, Only used for Kmeans. (Background = 0)\n            #if seg_slice.sum() > seg_slice.size*0.5:\n            #    seg_slice = 1 - seg_slice\n            \n            # Watershed\n            labels = watershed_img(seg_slice)\n            \n            # Writes labeled segmented object to return images                     \n            segmented_slices[i,j] = labels\n\n    return segmented_slices.astype(np.uint8)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from skimage.measure import regionprops\n\ndef find_left_ventricle(images):\n    \"\"\"\n    Returns numpy array (4d)\n    Finds left ventricle from labeled segmented images\n    \"\"\"\n    \n    num_slices, time, height, width = images.shape\n    segmented_slices = np.zeros((num_slices, time, height, width))\n    \n    all_labels = labeled_segmented_images(images)\n    \n    # Iterate over all slices and whole timeseries for images\n    for i in range(num_slices):\n        for j in range(time):\n            \n            labels = all_labels[i,j]\n            min_dist = 50\n            min_dist_label = 0\n            segment_found =  False\n            \n            # Iterate over every label in watershed labels to predict which is the left ventricle.\n            for label in np.unique(labels):\n        \n                # yx coordinates for labaled segmentation\n                yx_coord_labels = np.where(labels == label)\n                \n                # Do not count small or big segmatations (removes dots and background)\n                if len(yx_coord_labels[0]) > 8000 or len(yx_coord_labels[0]) < 100:\n                    continue\n                \n                # Upper right middle coordinates\n                cx = 3*(height/4)\n                cy = width/4\n                \n                # Calculates euclidiean distance between mean coordinates for segmentated labels and middle of image\n                euclidiean_dist = np.sqrt((int(cy)-np.mean(yx_coord_labels[0]))**2+(int(cx)-np.mean(yx_coord_labels[1]))**2)\n                \n                # Gets min distance\n                if euclidiean_dist < min_dist:\n                    \n                    # Check if segment shape is round.\n                    regions = regionprops((labels == label).astype(int))\n                    props = regions[0]\n                    y0, x0 = props.centroid\n                    orientation = props.orientation\n                    x1 = x0 + np.cos(orientation) * 0.5 * props.major_axis_length\n                    y1 = y0 - np.sin(orientation) * 0.5 * props.major_axis_length\n                    x2 = x0 - np.sin(orientation) * 0.5 * props.minor_axis_length\n                    y2 = y0 - np.cos(orientation) * 0.5 * props.minor_axis_length\n                \n                    d1_dist = np.sqrt(abs(x0-x1)**2+abs(y0-y1)**2)\n                    d2_dist = np.sqrt(abs(x0-x2)**2+abs(y0-y2)**2)\n                    \n                    # Checks if segment is round.\n                    if abs(d1_dist-d2_dist) > 20:\n                        continue\n                    \n                    min_dist_label = label\n                    min_dist = euclidiean_dist\n                    segment_found = True\n            \n            # Checks if we found a image or not\n            if segment_found:\n                # Writes segmented object to return images                     \n                segmented_slices[i,j] = (labels == min_dist_label).astype(int)\n            else:\n                segmented_slices[i,j] = np.zeros(labels.shape)\n                \n    return segmented_slices.astype(np.uint8), all_labels.astype(np.uint8)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Complete Preproc Pipeline\nHeavily inspired by the this paper: https://arxiv.org/pdf/1809.06247.pdf"},{"metadata":{"trusted":true},"cell_type":"code","source":"def preprocess_pipeline(patient, heart_pixel_size=150):\n    \"\"\"\n    [Patient Object] -> [4D np.array] (segmented left ventricle)\n    \n    Preprosessing pipeline for patient:\n        1. Rescale images (1 pixel = 1 mm)\n        2. Histogram Normalize (some images are brighter than others)\n        3. Crop images aroind ROI (identified using Fourier Transform over time)\n        4. Rotate images (such that left ventricle is in top right part of img)\n        5. Segment out left ventricle (for each 2d slice)\n    \"\"\"\n    \n    # Rescale images such that 1 pixel = 1 mm\n    rescaled_imgs = rescale_patient_4d_imgs(patient)\n    \n    # Histogram normalize\n    normalized_imgs = histogram_normalize_4d(rescaled_imgs)\n    \n    # Crop around ROI\n    cropped_imgs = crop_heart(normalized_imgs, heart_pixel_size=heart_pixel_size)\n   \n    # Rotate images\n    rotated_images = rotate_images_210_deg(cropped_imgs, patient.orientation)\n    \n    #return rotated_images\n    \n    # Segment out the left ventricle\n    segmented_left_ventricle_4d, labels = find_left_ventricle(rotated_images)\n    \n    return segmented_left_ventricle_4d","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Volume analysis"},{"metadata":{"trusted":true},"cell_type":"code","source":"def smooth(y, box_pts):\n    box = np.ones(box_pts)/box_pts\n    y_smooth = np.convolve(y, box, mode='same')\n    return y_smooth\n\ndef volume_for_patient(patient_images, slice_dist):\n    \"\"\"\n    Return numpy array\n    Array of total volume at each time for segmented images\n    \"\"\"\n    \n    num_slices, time, height, width = patient_images.shape\n    volume = np.zeros((time))\n    \n    if slice_dist == 0:\n        print(\"WARNING! Slice ditance is: 0 \\n Setting slice distance to 10.\")\n        slice_dist = 10\n    \n    for i in range(time):\n        time_volume = 0\n        for j in range(num_slices):\n            xy_size = np.sum(patient_images[j,i])\n            time_volume = time_volume + xy_size * slice_dist\n            \n        # Volume in ml instead of mm^3\n        volume[i] = time_volume/1000\n        \n    # Smoothing volume with convolution and removes last elements\n    #smooth_volume = smooth(volume,4)[2:-2]\n    return volume #smooth_volume","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Viz"},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\nimport matplotlib.pyplot as plt\nfrom skimage.util import montage as montage2d\n\nmontage3d = lambda x, **k: montage2d(np.stack([montage2d(y, **k) for y in x], 0))\n\ndef plot_patient_slices_3d(patient_slices, title=False, figsize=(20, 20)):\n    '''Plots a 2D image per slice in series (3D in total)'''\n    fig, ax = plt.subplots(1, 1, figsize=figsize)\n    image = montage2d(patient_slices)\n    if title: ax.set_title(title)\n    ax.imshow(image, cmap='bone')\n\n\ndef plot_patient_data_4d(patient_data, all_slices=False, num_slices=[0], figsize=(20, 20)):\n    '''Plots a 3D image per time step in patient data (4D in total)'''\n    if all_slices:\n        # Number of slices is equal to the first dimension of the patient image array\n        num_slices = range(patient_data.shape[0])\n    for i in num_slices:\n        plot_patient_slices_3d(patient_data[i],\n                               title=('Showing slice %i' % i))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Load ground truth"},{"metadata":{"trusted":true},"cell_type":"code","source":"import pandas as pd\nimport os\nimport numpy as np\n\nimport warnings\nwarnings.filterwarnings('ignore')\n\ndef export_patient_volumes(patient_ids=range(1,501), n_patients=500):\n    path = os.path.join(\"..\", \"input\", \"train.csv\")\n    ground_truth = pd.read_csv(path)\n\n    true_min = np.array(ground_truth.Systole.iloc[:n_patients])\n    true_max = np.array(ground_truth.Diastole.iloc[:n_patients])\n\n    min_vols = np.zeros(n_patients)\n    max_vols = np.zeros(n_patients)\n    patient_mask = np.zeros(n_patients)\n\n    print(\"Processing %i patients...\" % n_patients)\n    for pid in patient_ids:\n        i = pid - 1\n        try:\n            p = load_patient(pid)\n            if p != None:\n                lv = preprocess_pipeline(p)\n                v = volume_for_patient(lv, p.dist)\n                min_vols[i] = v.min()\n                max_vols[i] = v.max()\n                print(\"Sucessfully processed patient #%i\" % pid)\n                \n                # Mark patient as used\n                patient_mask[i] = 1\n    \n                 # Clean up data no longer in use\n                del p\n                del lv\n            else:\n                # Mark patient as unused\n                max_vols[i] = -1\n                min_vols[i] = -1\n                print(\"Error: Could not process patient #%i\" % pid)\n\n\n        except ValueError:\n            print(\"Error: Could not process patient #%i\" % pid)\n\n    print(\"Finished processing %i patients!\" % n_patients)\n    \n    print(\"Saving patient data...\")\n    np.savez(\"min_vols\", min_vols)\n    np.savez(\"max_vols\", max_vols)\n    np.savez(\"patient_mask\", patient_mask)\n    print(\"Patient data saved!\")\n    \n    return min_vols, max_vols, true_min, true_max, patient_mask","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Save results to output"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nids = range(1, 500)\nall_min_vols, all_max_vols, all_true_min, all_true_max, patient_mask = export_patient_volumes(ids)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Apply mask to filter out \"bad\" patients"},{"metadata":{"trusted":true},"cell_type":"code","source":"patient_indices = patient_mask.nonzero()[0]\nmin_vols = all_min_vols[patient_indices]\nmax_vols = all_max_vols[patient_indices]\ntrue_min = all_true_min[patient_indices]\ntrue_max = all_true_max[patient_indices]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Compare calculated min and max volumes to ground truth"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, (ax1, ax2) = plt.subplots(1, 2, figsize=[12,5])\n\nax1.plot(patient_indices, true_min, '*', label=\"true\")\nax1.plot(patient_indices, min_vols, '+', color='r', label=\"predicted\")\nax1.set_title('Predicted min volume for each patient', fontsize=16)\nax1.set_ylabel('Predicted min volume')\nax1.set_xlabel('Patient')\n\nax2.plot(patient_indices, true_max, '*', label=\"true\")\nax2.plot(patient_indices, max_vols, '+', color='r', label=\"predicted\")\nax2.set_title('Predicted max volume for each patient', fontsize=16)\nax2.set_ylabel('Predicted max volume')\nax2.set_xlabel('Patient')\n\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Correlation between ground truth and predicted values"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, (ax1, ax2) = plt.subplots(1, 2, figsize=[12,5])\n\nax1.plot(true_min, min_vols, '*')\nax1.set_title('Prediced min vs True min', fontsize=16)\nax1.set_ylabel('Predicted min Volume')\nax1.set_xlabel('True min Volume')\n\nax2.plot(true_max, max_vols, '*')\nax2.set_title('Prediced max vs True min', fontsize=16)\nax2.set_ylabel('Predicted max Volume')\nax2.set_xlabel('True max Volume')\n\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Correlation Coefficient\n# https://en.wikipedia.org/wiki/Correlation_coefficient\ncorr_max = np.corrcoef(min_vols, true_min)\ncorr_min = np.corrcoef(max_vols, true_max)\n\nprint(corr_max[0,1], \"Correlation predicted and true max values\")\nprint(corr_min[0,1], \"Correlation predicted and true min values\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Linear model for predicting volume\n* Least squares, 1 parameter - find scaling + bias to correct offsets\n* $y = ax + b$"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Create model, least squares fit\ndeg = 1\na_min, b_min = np.polyfit(min_vols, true_min, deg)\na_max, b_max = np.polyfit(max_vols, true_max, deg)\n\n# Use model\npred_min = a_min * min_vols + b_min\npred_max = a_max * max_vols + b_max","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Compare linear model to ground truth"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, (ax1, ax2) = plt.subplots(1, 2, figsize=[12,5])\n\nax1.plot(patient_indices, true_min, '*', label=\"true\")\nax1.plot(patient_indices, pred_min, '+', color='r', label=\"predicted\")\nax1.set_title('Linear model min volume for each patient', fontsize=16)\nax1.set_ylabel('Linear model min volume')\nax1.set_xlabel('Patient')\nax1.legend()\n\nax2.plot(patient_indices, true_max, '*', label=\"true\")\nax2.plot(patient_indices, pred_max, '+', color='r', label=\"predicted\")\nax2.set_title('Linear model max volume for each patient', fontsize=16)\nax2.set_ylabel('Linear model max volume')\nax2.set_xlabel('Patient')\nax2.legend()\n\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Ejection rate calculation\n* $ER = \\frac{V_D - V_S}{V_D}$"},{"metadata":{"trusted":true},"cell_type":"code","source":"def ejection_rate(vd, vs):\n    \"Returns a float between 0 and 1\"\n    return (vd - vs) / vd\n\ntrue_ejection = ejection_rate(true_max, true_min)\npred_ejection = ejection_rate(pred_max, pred_min)\n\nfig, ax = plt.subplots(1, 1, figsize=[12,5])\n\nax.plot(patient_indices, true_ejection, '*', label=\"true\")\nax.plot(patient_indices, pred_ejection, '+', color='r', label=\"pred\")\nax.set_title('Ejection rate, true vs pred', fontsize=16)\nax.set_ylabel('Ejection rate')\nax.set_xlabel('Patient')\nax.legend()\n\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Calculating MSE\nmse_error = sum((true_ejection-pred_ejection)**2)/nbr\nprint(mse_error, \"MSE\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"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.6.4"}},"nbformat":4,"nbformat_minor":1}