{"cells":[{"metadata":{"trusted":true},"cell_type":"markdown","source":"# Images for presentation and exam"},{"metadata":{"trusted":true},"cell_type":"code","source":"import pydicom as dicom\nimport re\nimport os\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom skimage.util import montage as montage2d\n\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":"### Load Patient(s)\n\n"},{"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":"### Plotting functions"},{"metadata":{"trusted":true},"cell_type":"code","source":"def 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, slices=[0], figsize=(20,20)):\n    '''Plots a 3D image per time step in patient data (4D in total)'''\n    if all_slices: \n        slices = range(patient_data.images.shape[0])\n    for i in slices: \n        plot_patient_slices_3d(patient_data.images[i], \n                               title=('Showing slice %i' % i))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Fourier Transform"},{"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 get first harmonic mean\n    fft_img_2d = np.fft.fftn(image_3d)[1, :, :]\n    return np.abs(np.fft.ifftn(fft_img_2d))\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    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, threshold is used for finding a fitting threshhold for segmentation\n* We also tried k-means but with less good result\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\nimport pandas as pd\n\nfrom skimage.morphology import binary_dilation, binary_erosion, binary_opening, binary_closing, disk\nfrom skimage.filters import threshold_otsu\nfrom sklearn.cluster import KMeans\n\n\ndef kmeans_segmentation(patient_img):\n    #Code for kmeans\n    xx, yy = np.meshgrid(np.arange(patient_img.shape[1]),np.arange(patient_img.shape[0]))\n    patient_df = pd.DataFrame(dict(x=xx.ravel(),y=yy.ravel(),intensity=patient_img.ravel()))\n\n    km = KMeans(n_clusters=2, random_state=2018)\n\n    scale_patient_df = patient_df.copy()\n    scale_patient_df.x = scale_patient_df.x/250\n    scale_patient_df.y = scale_patient_df.y/250\n    scale_patient_df['group'] = km.fit_predict(scale_patient_df[['x', 'y', 'intensity']].values)\n    seg_pat_img = scale_patient_df['group'].values.reshape(patient_img.shape)\n\n    return seg_pat_img\n    \ndef thresh_segmentation(patient_img):\n    \"\"\"Returns matrix\n    Segmententation of patient_img with threshold\n    \"\"\"\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 segment_multiple_kmeans(patient_img):\n    \"\"\"Returns list\n    List of segmented slices with function kmeans_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 = kmeans_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    \n    seg_slices = segment_multiple(patient_img)\n    num_slices = patient_img.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\nEqualize_adapthist: (Gave better result) \nAn algorithm for local contrast enhancement, that uses histograms computed over different tile regions of the image. Local details can therefore be enhanced even in regions that are darker or lighter than most of the image. Region = 1/8 * image size\n\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":{},"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\n\nfrom skimage.measure import label\n\ndef labeled_segmented_images(images, kmeans=False):\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            if kmeans:\n                seg_slice = kmeans_segmentation(images[i,j]).astype(bool)\n            else:\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)\n\nfrom skimage.measure import regionprops\n\ndef find_left_ventricle(images, kmeans=False):\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, kmeans)\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 = 75\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 upper right corner 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                    # This should be d1_dist/d2_dist instead...\n                    if abs(d1_dist-d2_dist) > 30:\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)\n","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":"# Example patients\n\nWith patients 100, 150, 400"},{"metadata":{"trusted":true},"cell_type":"code","source":"patient_100 = load_patient(200)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 1. Rescale images \nRescale images to 1 pixel = 1 mm"},{"metadata":{"trusted":true},"cell_type":"code","source":"rescaled_patient_100 = rescale_patient_4d_imgs(patient_100)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Plot all time steps for slice 5"},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_patient_slices_3d(rescaled_patient_100[5])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Plot all slices for time step 2"},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_patient_slices_3d(rescaled_patient_100[:,2])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 2. Histogram Normalize\nSince some images are brighter than other and we want same contrast over all images we do histogram normalization.\nHere it would be good to check intensity normalization also?\n\nImages seemed to not be noisy and no noise filter was used.\n![image.png](attachment:image.png)**","attachments":{"image.png":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAlgAAAENCAYAAADTzRsFAAAgAElEQVR4Ae3dCbQ0Z13n8W8nb97skJAAISGsgbAm7HsMmywCggJDBFREwI0ZBgV0QAWPIszgqMgqAzOMAsPAHAHZZGSQxYDsm4R9URZJAiSEkOTN8j5zft1V6bqd7nv73lt9b3fV9zmnT1V3Vz9d9am6df/9rHDVdGPgfcCvA4Orvu0rCiiggALAXuANwF8CByiigAIKNAWm3RROBvK4G7CnubHrCiiggAJXCtwIuClwO+BqV77qigIKKDDlV1cCrlsB1wJuCxyjkgIKKKDAVIFbA9cHjgeyblJAAQWuFJgswToauGNVNXht4A5XbumKAgoooEAtcCBwKnB14BrAXeo3XCqggAIRmAywjqsCrLyXIu9UE5oUUEABBdYK5AdoqgbTTvWQav2wtZv4TAEF+izQDLByo7glcEIFclD1Cy2/zkwKKKCAAmOB3CcTYNXpZkA6CJkUUECBoUAzwEpAdS8gRd91ugmQG4dJAQUUUGAkkPtmqgfTVrVO3itrCZcKKDAUaAZYR1ZVgs2hGdKAM6VaJgUUUECBkUCqBO85MYzN4VWJ1sEiKaCAAhFoBljpBVNXD9Y6GeclvQlz8zApoIACCsBRwJ2nQNwVSEchkwIKKLAmwDoNOGKKSW4kzaLwKZv4kgIKKNAbgfzoTCP3yXSbKT9SJ7fxuQIK9ESgLsFK75cEUtOKt28O3LAnHh6mAgoosJFA2qoeOmWj/EC9+0TNwJTNfEkBBfogUAdYacg+K4hK0JWbhqO69+GK8BgVUGA9gYx7dfsZ98PcTxN8ea9cT9D3FOiJQB1gnQJcb8YxZ5s06JxWujXjI76sgAIKdFKgHr192sGlg1Dez8juJgUU6LlAgqcETmk7MK39Vc2TEq4b1E9cKqCAAj0VSPurDMg8Kx0LpLG7SQEFei6QAOu6VU/BiiIl4I9b+9JoOohUE5oUUECBvgqkrWp+jFbtr1IT+LPA/ZsjNuSHagKs5nA3ffXyuBXovcB9gPOAAkcUeFaB8wq8r8CpBQYF2A+82sabvb9WBFCgzwIpyf/Y6F55YIFHF/hmgS8WeEiBvJb7KO8BrtNnKI9dAQVGAr9X3TAuht+6AH5UoFSP/7cPbvrD6qbx4XUawmupgAIKdF3gp4CLYHA5PPgc+FbjXvmFAvfOD9UrgG8C+eFqUkCBngt8APghHPNC+MGnx8HVMMi6EF74MuAjwDnAQ3tu5eEroEA/BTKV2HOAy+CwN8JZr4Oyf+398j0vh8FbqhL/p1hN2M8LxaNWoCmQwOl58PjjoZy59oZRLoLyKOAewGeBlHbVPQ+bebiugAIKdFngGOD/Aq+B250E5S8m7pUpzXo67M24ge8E/orRiO9dNvHYFFBgHYEESy8aBVivSvH2tOApjTX/Efit6hdZbjQmBRRQoE8CGbn9LOB34OP/OqN0agCXfr66V57tDBh9ujw8VgWuKpBuMM8fFXtvON9gfr19Dbj8qtn4igIKKNBZgfzIvBB4LnDu+mMClgEMEoi9HDi/syIemAIKbCiQAOuyDbcab/CV8aprCiigQC8E0jMwpVZ1mmcIhq/WG7tUQIF+CkyrEuynhEetgAIKKKCAAgq0JGCA1RKk2SiggAIKKKCAArWAAVYt4VIBBRRQQAEFFGhJwACrJUizUUABBRRQQAEFagEDrFrCpQIKKKCAAgoo0JKAAVZLkGajgAIKKKCAAgrUAgZYtYRLBRRQQAEFFFCgJQEDrJYgzUYBBRRQQAEFFKgFDLBqCZcKKKCAAgoooEBLAgZYLUGajQIKKKCAAgooUAsYYNUSLhVQQAEFFFBAgZYEDLBagjQbBRRQQAEFFFCgFjDAqiVcKqCAAgoooIACLQkYYLUEaTYKKKCAAgoooEAtYIBVS7hUQAEFFFBAAQVaEjDAagnSbBRQQAEFFFBAgVrAAKuWcKmAAgoooIACCrQkYIDVEqTZKKCAAgoooIACtYABVi3hUgEFFFBAAQUUaEnAAKslSLNRQAEFFFBAAQVqAQOsWsKlAgoooIACCijQkoABVkuQZqOAAgoooIACCtQCBli1hEsFFFBAAQUUUKAlAQOsliDNRgEFFFBAAQUUqAUMsGoJlwoooIACCiigQEsCBlgtQZqNAgoooIACCihQCxhg1RIuFVBAAQUUUECBlgQMsFqCNBsFFFBAAQUUUKAWMMCqJVwqoIACCiiggAItCRhgtQRpNgoooIACCiigQC1ggFVLuFRAAQUUUEABBVoSMMBqCdJsFFBAAQUUUECBWsAAq5ZwqYACCiiggAIKtCRggNUSpNkooIACCiiggAK1gAFWLeFSAQUUUEABBRRoScAAqyVIs1FAAQUUUEABBWoBA6xawqUCCiiggAIKKNCSgAFWS5Bmo4ACCiiggAIK1AIGWLWESwUUUEABBRRQoCUBA6yWIM1GAQUUUEABBRSoBQywagmXCiiggAIKKKBASwIGWC1Bmo0CCiiggAIKKFALGGDVEi4VUEABBRRQQIGWBAywWoI0GwUUUEABBRRQoBYwwKolXCqggAIKKKCAAi0JGGC1BGk2CiiggAIKKKBALWCAVUu4VEABBRRQQAEFWhIwwGoJ0mwUUEABBRRQQIFawACrlnCpgAIKKKCAAgq0JGCA1RKk2SiggAIKKKCAArWAAVYt4VIBBRRQQAEFFGhJwACrJUizUUABBRRQQAEFagEDrFrCpQIKKKCAAgoo0JKAAVZLkGajgAIKKKCAAgrUAgZYtYRLBRRQQAEFFFCgJQEDrJYgzUYBBRRQQAEFFKgFDLBqCZcKKKCAAgoooEBLAgZYLUGajQIKKKCAAgooUAsYYNUSLhVQQAEFFFBAgZYEDLBagjQbBRRQQAEFFFCgFjDAqiVcKqCAAgoooIACLQkYYLUEaTYKKKCAAgoooEAtYIBVS7hUQAEFFFBAAQVaEjDAagnSbBRQQAEFFFBAgVrAAKuWcKmAAgoooIACCrQkYIDVEqTZKKCAAgoooIACtYABVi3hUgEFFFBAAQUUaEnAAKslSLNRQAEFFFBAAQVqAQOsWsKlAgoooIACCijQkoABVkuQZqOAAgoooIACCtQCe+oVlwoooMD8AmUA5P5xLeA44AjgEODQxmPyeba/DLi48bhkyvqPgH8DzoXBFfPvk1sqoIACyyNggLU858I9UWDJBMqBwMHA3ipoui5wI+DG1eN6wDWAI6ttcj85qHo015PPZLq8CrbqZQKven0fkCDr+1C+AXwN+Gr1SOCVoOxSYB8M9k9m7HMFFFBgGQQMsJbhLLgPCiyNQEmwlBKpBFO3AW4LnAQkmDq8EUAl6NpOyr1nnvtPqQKxBFQJwhJ4/QvwZeDjUD4LfHtU4jVIyZhJAQUUWAqBeW5wS7Gj7oQCCixCoKQd5onAzavHrYFbVeup9tvtlKrIBHN1QHd0FeydBjweOB/4HPDPULI8q3p8FwYJzkwKKKDArggYYO0Ku1+qwG4KlFT7pYTqvsA9gJtUz48HEtC0nRLo1I/kXz/a+J6jgLtXj1Qxfgf4JvAlKP8A5HEODFICZlJAAQV2TMAAa8eo/SIFdlOgpG3UCcApwMOA04FrV9V+W92xtJW6sKq2S/VcApy6HVWe/7h6ZJusJ8hJw/dUNTYfaRife1H9yPOUnqW6MsHgvCmfT1VmHgm6HlFVH74bylurkq3v2G5rXk63U0CB7QgYYG1Hz88qsPQCJb38bl+VVN0HOLUKWjZTUpWG5N8Dzh6VBg2X362Cl3+tSozOqYKoBFYJvJq9/6ZV1TW/P43g6x6HCbyuWVVbJlBKSVvahOWRgDDHc+yc7AnQbgacXFUnfgR4D5QPjtpvDVK9aFJAAQUWImCAtRBWM1VgNwWG7apS3ZcSnAdUpVbXmXOPEkylFOo84JPAx4DPV4HV90c9+xJsDS6aM795NkswltKtC6qN02vww+MPlsOqoCuB1TFVoJW2Yneoju1qVelXM2gbf3xUJZk87lk9EhR+EsrbgLeMgkfbazXBXFdAge0LGGBt39AcFFgSgZK/5xsCPw08tiq5SXXbRinVdwmovgWklOdM4DNVMHU+DNJ7bxfTMJhLz8E8qjRsR5b2Vyntul1VJZiSugSSaQi/3nHX1Yj3rkq2/ieUd46O32EfamGXCiiwPQEDrO35+WkFlkBgGFil59+DgUdVPQCnjT3V3NeUGn2lKqH6R+CjwBeqsaVSgrXkaZBqyFRZnl31Hnxd1a7rFsBdgLtVgVfG7ZqVUoV416oK9VPA66G8CwbpiWhSQAEFtiVggLUtvp37cIF0p88/hLShyUCP6baeQRc/nbYwg1G1zs7tkN+0BALDqsA0Wn80cP+qumy9/Ur13w+BD4zaIvFPo+tnkIE7VzgNq/cSFObYPjR6lFdW43gl0ErbswRSad81rRoxf0t3qgKtn6uqDl8/GmvLqsMVvjDcdQV2VcAAa1f5N/7yAmk7kqqPh1YNldMTLG1OEnClHUwaH3+mQNqT/P1g2CV943zdYtUFShp/P6aq4spAoOvNK5q2TRkJ/e1AqsK+AoM0Su9wGmRA0g9ASencX1WN3R8C3A+4QRVsTR5/Sv3uWAVmPwO8AsobYZC/MZMCCiiwVYFyOJQPQSmNx0VQzthqjn5u6wIF9hS4Y4EXFfhugcubJ2bK+o8LvKvAw8to+pKtf7mfXGKBcjSUlLIkeLi08bc65ZIo50J5C5RHQzkWhlWJS3xsi961DFVRjofyRCh/D+X8DfwuhvJ2KA+B0hh0tRwC5UVTPvsMGM7RuOgDMX8FFFgtAQOsZTlfBa5X4DkFPj/tv+YGr51X4LUFTq+qFZflsNyPbQkMg4PTofz1HIFB5vB7HZSHQ8lcgW2lVK/NerT1HXU+s75nWhVf/ZlNLMu1ofw8lL+BcsGUYKn5Z/ZdKC+DcgdItWwa2BtgbQLbTRXopYBVhEt22gvcFHh+1WA5g0NuNqVnVdrkpAv7swq82fZZmyVctu1LqoV/vTqvqd6alc4F/g547ajR+uAHszbcwusZS+r3q955qXLMeFf1/IBp/1Q/6kmb87y5Xr+fr859Z9Yj1XSz3rt61fbwGUAmgd5GGqSBfILVd1WN4h9XDb46LSDN+Fu/Wg3x8D9GjeGHI9Nv4/v9qAIK9EjAEqzdPtkFblfgPQX2N38+b2P9OwWeUNbvsr7bh+33zxQYllrdG0pGIl+vOjAlVv8LSkq40mZvESkjv6cxfD3lTRrMpydiHTjVAVWCqjwSfKWnXx4JxtJesB7ZPevJq34/29bBWv35ZnCW78kj350xudbrGbjFY88k1+UBVZXqD9cp0bpk1Ag+Y2iV/RPbWUW4RX0/pkAXBSzBWoKzWlXl5R/Y84A7t7hLGRPouRkXqMArBqNeVi1mb1aLExiOwP5LwK9UY1tN+6oENe8HUqrytzCoB+qctu12X0sbpHrC5eRVV+El6EmwlACpDrISfNWBWLbNAKHNcanyXkZRT8CVxvn1I6VX+Y76Ma068LPVZ7d7PBOfHzaK/zsoGQcsA7SmRCs9CyeHu8jUPT9VHe+0/ZvI16cKKNBXAQOsXT7zZfSPKqNt/xfglgvYnUwt8qwqyHr+YDxa9gK+yiy3LzAceiEDZ/4O8EBGvUgns00Ak9HOMxTB/4bBNqvLJrOf+jxBUh1Q5PtePRrGYDg0Qkqi6lKmLOsAKxklMMu1fdtGrgkEUw2e4LAO1LKsqwfrQCsDhp4IpEdfqryzzZeruQ8b2bW5OqxWTe/B91aDtSbITRVtfez5sqxvZo7ENnfQvBRQYPUErCLcjXNW4C5bbMy+prvnHNWIFxb47bK2FGI3DtnvnCmQ6r3yS1A+P1H11Dy934OSAOA2O9grMAFP2j2l5CnDH/w8MG/7wARIKRWqS7SyTIA274+J/Aj8uWqk+QRuT5jJ1/obZS+Uu0B5LZT1qg3r8/O01nfBDBVQoAsCBlg7fRYLnFTg3fXdeQeWGe7hUfYu3OkzPc/3pUqw/BGUBFDTLoW09/kolMdAhgnY0ZTv+7MqSEqp07yTLWcnM65URohvBlgfrwbNnfcgblMNqptqxUwDtMOpXA3Kk6D884xzU5+vV7fca3OHj9OvU0CBBQkYYC0Idmq2BY5Ju6gCl9Z35x1afra0285r6vH54mYEyslQMh9exl2adhlcCBmZvGSC491IqebLyOYJklLdvJmmBWmv9O1GgJVSqDduMo+bVHlksNSMzL5LaVia9YZ1ztOPoLwEhr0+d2kf/VoFFFhCAQOsnTopaXdV4IkFMjjo1P+oC3799WVzpRA7RdOz7xmOqZSef++FcsWUSyGlVl+sBsbM6P27lTL0R6bXyVQ0993kTqSx+IWNACsN4f9gg5HnJ78iQ0RkWqhUNWZ9F1M5Bkp6C84qabys6omYdnQmBRRQIAIGWDt1HRS4VYGzFhxErRe4XVzgV6sG9jt12H7PGoHhEAwZkf1LUwKrXBoZliGDYN5pCUYHvybw9Wpi6JutOYyNn6TtVt3oPSVgaRCfCambjcY3yiUld5naJ2N8pdPGLqfhvfI1M85b/Wf9MSgPWoJzt8tWfr0C/RVYb/6y/qos8MjLqFfVf2Q0oOgCv2ndrNOm5slzTA68bia+uVWBBFfDhtsvAFL9NZkyfMGLgafA4COjAs/JTXb0eSZJTgnaWVVJ0rxfnqrEDNLZDKYynMOXqhKtefNJI/v0TkyQdd68H1rgdimF22g/Mn9o2q09zCBrgWfCrBVYDQFLsBZ9nqqqwTQy/179M3cXl1dUbcDyz9O0YwKZD3DYU/DbM0pA8vpTIANfLk3KGFb3rILBZrC00Q6majGjyjcbuCe4Ss/CzaQMnvoTu189WO/yzLkIp/05fxnKI6AkSDQpoEA/BQywFn3eCxxX4K3T7sK79Nq3Czxo0cdt/rXAsNv/L0M5e0Zw9UkoP7ODwy/UO7ao5fWB90wEWG8DjlnUF+5MvjMDrH1TRnfPn/ZXoTzKIGtnzo7fosCyCFhFuENnohoa4d7VfGc79K0bfs3xwL8rkJIG00IFMkEwv1iNrD+tHVEGtvx1GLwJBqmC6kJKteJxEweSgUIzVU4XU4LHj005sEztk4FVU5LlPXcKkC8p0EUB/9h37qxmEtnHb3L8n53Yu4wr1Ob0PDuxzyv2HSm5Gk698kdVm6Tm/qf67N3AU2HwoeYbHVjP5MyZrqmZvlJNrdN8rSvr/wSkfWV6O06mTNKdIOsMS7ImaXyuQDcFDLB27rxmPKBdHMNn5oGm9OpJBWyLNZNoO28MS64y+njmhJwsuUpw9S7gN2Hwqe18y5J+NsFVgqw6peTqXzfZwL3+7CosBzD44KhzAlnm/DZTgqzMN5rBYjczllgzD9cVUGBFBAywduBElVEPrEdPTHi7A98891ecVs31NvcH3HAegZLG4QmuMu7TtHZHGXbgaTDIBMZdSwkgUjXWbBT/naonYNeOtXE8JUFWSrJ+E6YGWdergu1HG2Q12FxVoIMCBlg7c1JTcnXqznzVlr4l1ZcPdwqdLdnN+NCw5OqJVXA1ObVMSjbeUQVXn5uRwaq/nGrRG08cREZ0z1ALPUiDD4+qfUm172RJ1nWrICslWfYu7MHV4CH2U8AAa8HnvZpcOT31Mh7Qsqbc5NMN/hbLuoOrtV/DhswZTPP3p5RcZdDNvwV+CwYZV6qrKQHWSRMH16MAK0c++GhVXZhR8KcFWSnZfLAN3yeuEp8q0BEBA6zFn8ibVz0Hm1Uli//WzX/DTYF7Obr75uGmfCJjRk0LrvJP9i3A02GQCZC7nFI9mirCOuXYvwVkENUepUF6Faa6ML1EJ4OsDGPxbMBpdXp0RXio/REwwFrgua6q3O4BbHZ6kQXu1cys8w/xflMaYs/8gG9MEyi3BJ4zpXosJVdvAp4BgwxV0PWUBu7N4T8uBr7R9YOefnyDj1dB1vunBFm3HV0vJdWGJgUU6JCAAdZiT2Z6UCVoydQoq5DSVmza1C2rsO9LsI8lQcXvAnef2JkEVxkjKcFVhinoQ0qJaPO6Px/4Wh8OfPoxDnuJPrVq+D65yQOAZ0Jp9ric3MbnCiiwYgIGWIs9YZkS5C6L/YpWc09j93sUsAv5pllLpnNJVdDDgcm/q5RgPAcGX910tqv5gVSHTwZYF/S3BKs+iYNPVlWCkyWYaQP5WOBXYdjztP6ASwUUWGGByX8EK3woy7XrVVumNBxP0LJKKSVuCRZMcwsMe4L9UsYTmyi1SQ7fHPUkHP5znTvHDmx4MmsD9ZRgpQ1W31PaYmXA2e9PQGTuyacAD5nS6D0BmPfqCTCfKrDsAv7RLu4MpRdViv5XrTTo1lXpw+JkOpVzxj0ig8g+ndF4Z82j+yHwx0DGu+pTyqC1GVSz2bEjA4ymHVbP0+AK4A3AC6eMaJ8q5rTfu+ME0ikTHQYm3vapAgoso4AB1uLOSn7B32px2S8s58wf95MLy717GaeX6DOB9AhrpgQTLwFeA8N/qs33ur6ewTSPbhxk5lZMtdhkL7rGJn1aHVwCvAx4LXDZxJHneno2lBMar6ejzH0txWqIuKrACggYYC3uJJ0OXHNx2S8s55S8/YRT58zjW1Kt8xvAnSa2TinFm4EXweDCiff68PSGEz0IE2Cl/ZkB1pVnf/C9am7CVBlOptw7fhnKIcARwO2rayw9fU0KKLAiAgZYCzhRVSPx3BRXtS1T/kGmkbJppsCwavAho3nlrtI+JpP9/gEMvjvz491+Y1qA1Zfek5s4s8Mepb8HfH7iQ7lvpD1fevWmo0yqDDP8hwHWBJRPFVhmAQOsxZyd3BQnR7FezDctJte0BUlbLNNsgVTl/PbEZMbZOg2502Pwi7M/2ul3ck9JgNUMBupJnjt94Fs8uATjfwikE0AznQDlmXC3VA3mx06mHVq1DjPN43Fdgd4JGGAt5pQnuMo/mVVNaYd1iwLOkzb1DJa0L8qYRpNTC6Xd1X8DMjVKX1PGcsoPjGZKA/eejeDePPyp66mKPxYG6Qjw1lFbPTJeWiPtOw2u9eSqo0xKtSavt8a2riqgwLIJGGAt5oxksM5lnntwnqPO6PPHzLNhv7YZzjP4QOARU3qIJrB6FQz63FvuWkAauTdT2l/ta77g+rA92n1GvQkHD4KXvh0u/CSURju1C/bCl+qq+gzaepuJnpkyKqDAEgsYYLV8csqoUWq6Va966U9+LTt9x1Wvj3Rc+IWJRtzZKlWD/xkGmdC4z2lagJUR3A2w1l4V6T34D0A6RLwUnvxSeMAR8KrL4NNA2sBnrugMozZMGe4lAdaq31fq43GpQOcFVm2MplU4IWknkRvhqqcMO5CxjD6x6gfS8v6ni/3fA8dXw3Ckiif/LF/e86rBmjnDCyTIaqYEWJPDETTf7+v6OcBzR+3VyuPhzD3woerSygQQP5gcOizNDlKqfHZfwTxuBVZJwBKs9s9WqgbTAHrVU9qInFLWzie36sfUwv4PMnjon1WlWC8GzgMyie9fw6DvQUSqsTL+W7OU5VIgpXqNqq8WTkN3sji3mj7nNaPSrDTDSmHo/wHeA2SEiytT2v7F16SAAisgYAlWiyepmh4n7a8yknUXUgZKTePaBBWmKwUG+S/4KSi/UzVQTgPuK+tyrtysfyu5VjKcQDOllCZFMabZAhnOI8M15H58xpS2ffUnE2ClbWQCepMCCiy5gAFWuycov9xTetWcIqTdb9jZ3PJr+eCd/cpV+rZBhh9IdaFpJJAfFpMBVoIHA/SNr5AUW/1uNZflI2eM2p5BR/M3mfuLJYIbm7qFArsqYBVhu/wJsNI4vCuuaU9zbLtE5tZhgVSPT3aMSHshA6z5Tvq/AM+qSkVnBVBph5WhMEwKKLDkAl0JBJaFOaU9dbfqZdmn7exHpuro0vFsx8LPri+QUpUMTptrpplSRXhB8wXX1xXIkBapen4HTI6LNfxcOp9MBrHrZuibCiiwOwIGWO26Z4DFLo22nCrkDDjalSrPds+2uTUFco3cYUr7oYw3kKpU0/wCX6hmCUj182RJVgKs5kTQ8+fqlgoosKMCBljtcqe0Jw19u5ISYKVRrQFWV87o4o4j7a9uN6V6PNWDk0HC4vaiOzl/rirJmmzQnh9wGT7FpIACSy5ggNXuCUoD1OYcbO3mvvO5JbDKtD+T1T47vyd+4zILpGo8c+alB+1kOhLIkB+mzQt8Cng6cGYjSM3fZDrSdOk+s3kZP6HACgjYi7Ddk5R/MF3rdZeBDVP12dfJi9u9QlY/t/yDTzVVpsNJqVU6QWQ4j58FMsr9ZHpUNZjTZ6rG7hnRPVWGX5oywfHkZ30+Gs49QVbGXrtzBRLvBK59npLJa0OBpRcwwGrpFJVR1eB1WspumbLJP9E0qjXAWqazsnv7ktKopwGPATIp+Eal4KnOSs+4pEwLk+AqY4b9BvDe6nUX6wtkePdnAH8O3LYqwUpPwnQgMCmgwJIKGGC1d2JS0pN/OF1LaVN2XNcOyuPZskDaU325Gmo8wVIemT4oQ47nkSAqj5R05f6S0d0TlKVkN1XNuZ7ymb7P2bjZE5C2WCnJemHVszfDNeQ8mBRQYEkFDLDaOzGpKuni+DQpwTLAau86WfWcMvVN5l3MmG8Z0T6PBF11Q/Z6PQFWUpZ5pKSrXmYbJ3+ugDaxyNw5KT18AXB74N0zhnLYRJZuqoACixIwwGpPNiVYXQywUvqQASRNCtQCBke1xM4uE5i+qyoVzGzQuX8n4DUpoMASChhgtXdSulpFGKFjM+nzAPo+mXF7V4s5KbA1gQRZbwe+VlXFbi0XP6WAAgsX2KiB6sJ3oENfkCrC9OzpYupy8NjF8+UxdVsg1bIZJytt3UwKKLCkAgZYLZyYaqTzdFHvqmcGNzyqBSqzUEABBRRQoBcCXQ0IdvrkZdC/a+30l+7g9xlg7SC2X6WAAgoosPoCBljtnMP0tOt6gNXFBvztnP/aQ+QAAA+RSURBVH1zUUABBRRQYELAAGsCZItPE2B1uafd0VYRbvHK8GMKKKCAAr0UMMBq57R3vQQrg0SmEb9JAQUUUEABBeYQMMCaA2mOTbpeghWC4zNUwxwWbqKAAgoooEDvBQyw2rkE0gi8q0M01ELHd3Ai6/rYXCqggAIKKNCqgAHWNjmrIRpOqKYB2WZuS/3xTGSdqkKTAgoooIACCmwgYIC1AdAcb2dOtgRYXU8JsDJZr0kBBRRQQAEFNhAwwNoAaI63M93QiT0owbKKcI6LwU0UUEABBRSIgAHW9q+DBFjX3X42S59DRnJ3NPelP03uoAIKKKDAMggYYG3/LPQlwEpVaB8Cye1fEeaggAIKKNB7AQOs7V8CGR+qy6O410K5Vm5aLPWsPVwqoIACCigwU8AAaybN3G+c3JPG3ynBurnVynNfF26ogAIKKNBjAQOs7Z/8m/Vk+IIBkGNNlahJAQUUUEABBdYRMMBaB2fOtxJ07J1z21XfLPMtHrfqB+H+K6CAAgoosGgBA6xtCBc4FLhRD4ZoqJWOAG5SP3GpgAIKKKCAAtMFDLCmu8z7asa/OmbejTuwXQKskzpwHB6CAgoooIACCxUwwNoe7/WBzEPYl3SYAVZfTrXHqYACCiiwHQEDrO3owQ17VoKVhu43KJBAy6SAAgoooIACMwQMsGbAbPRygYOqXnV9CzZSaneDjXx8XwEFFFBAgT4LGGBt/exncNGMC9W3lODqxn07aI9XAQUUUECBzQgYYG1Ga+22mfz4Fmtf6sWztDk72RHde3GuPUgFFFBAgS0KGGBtAa5A2iKlN90JW/j4qn8kx34KcLVVPxD3XwEFFFBAgUUJGGBtTTYDi94eyPQxfUyn9qz3ZB/Psce8eYH89iqb/5ifUECBLgo47cnWzurhwJ239tFOfCqDjV4X+FonjsaDUKAVgW9eDQ49ES7e9SDrGCjvhAPuOCptb+XozKQfAu+E/Y+FwQ9W/9pJbctFwPeB/btx9gywtqZ+vZ42cK+1DgbuXuDMAVxRv+hSgX4L/PcHwcV3ZdTDeDeDrMEAfnw55O/0kH6fE49+swL74JxLIYUIeaxySg3dmcCzgR/vxoEYYG1N/T7AkVv7aCc+lQv3fsCLgAs7cUQehALbFjgnP7xuvQwTohe44IBRoJfpvEwKzC2wH76xH44Grj73h5Z3w/N3symPbbA2eWGU0S/C+/ZogudZQhmioo+9KGd5+LoCKbXalaqIKfTZj2XZlym750tLLNCla2dX/wYMsDZ/lacH3cmb/1jnPnEUkEDTpIACCiiggAITAgZYEyBzPL07cJ05tuv6JmnfcbficA1dP88enwIKKKDAFgQMsDaBViClNuk9aMPRkVuqCVOiZ1JAAQUUUECBhoABVgNjjtU0YE0vIdNIINPmnF7AzhJeEQoooIACCjQEDLAaGOutFsjgovcG0lPINBLI9fOAakwsTRRQQAEFFFCgEjDAmv9SuDbwM/Nv3pst7wDcqZo+qDcH7YEqoIACCiiwnoAB1no6a9+7F3CztS/5rGqP9ijgMDUUUEABBRRQYCRggDXHlVBGA66d4dhXM7F+ArjdzHd9QwEFFFBAgZ4JGGDNd8Izanl6D2ZuI9NVBa4B/HI1COtV3/UVBRRQQAEFeiZggLXBCS9wTeDnGU0dsMHWvX0711EGHb1nbwU8cAUUUEABBRoCBlgNjMnVquF2Sq9Os/RqUucqz08AzqiqU6/ypi8ooIACCijQJwEDrPXPdkZs/0VGA4yuv6XvRuBBGRdLCgUUUEABBfouYIA14wqoBs/8BQOGGUDTXz4WeGqB46e/7asKKKCAAgr0Q8AAa/Z5TqP2X7Pn4GygGe9krsZfcXT3GTq+rIACCijQCwEDrCmnuUDaEz3NEcqn4Gz80kHA44Gf3HhTt1BAAQUUUKCbAgZYE+e1wBHAfwAeCOgz4TPn0wSov1vg1Dm3dzMFFFBAAQU6JWAA0TidjXZXqRo8uPGWq5sTyHhhmRQ7QZbtsTZn59YKKKCAAh0QMMCqTmIZBVTpMfgs4MgOnNvdPoQEWQ8F/qDAibu9M36/AgoooIACOylggAUUOBx4MvA8LHFp8/pLe6zHAX9SnMexTVfzUkABBRRYcoE9S75/C9+9Apnm5TerdleWXLUvnmvsEcA1CvynAXys/a8wRwUUUEABBZZLoLcBVtXe6k7AE4CfAw5ZrlPTqb1JSWmm0jmqwAuBtw/gvE4doQejgAIKKKBAQ2COAOuKtKVZ+UmO3wAHPHJ0HKm2un4mJwYeBpzU8HB1sQJ3AP4i7gVeCZwJXPJG2H8WlOcMa2sXuwMdyr106Fg8FAUUUKBzAhsEWJfvhxdn4Mj0qDtwVY/+MNjzX+HIY+HQe8EpQAYRzTQ4Gxz/qh7xUu/30cDDgXsBn/sRfPBtcN474SJGj6Xe+SXZuS9UwemS7I67oYACCigwKbBBgHHZfnjvY6sRzSc/uzLP85/7qzA4f7THK18atzLw6+9o2r6ddinc4+vAuZZera+19t1XGGCtBfGZAgoosGwCGwRY2d2SgMTehst25rqzP3XAWy+7c2SLOxKtFmdrzgoooEArAgZOrTCaiQIKKKCAAgooMBYwwBpbuKaAAgoooIACCrQiYIDVCqOZKKCAAgoooIACYwEDrLGFawoooIACCiigQCsCBlitMJqJAgoooIACCigwFjDAGlu4poACCiiggAIKtCJggNUKo5kooIACCiiggAJjAQOssYVrCiiggAIKKKBAKwIGWK0wmokCCiiggAIKKDAWMMAaW7imgAIKKKCAAgq0ImCA1QqjmSiggAIKKKCAAmMBA6yxhWsKKKCAAgoooEArAgZYrTCaiQIKKKCAAgooMBYwwBpbuKaAAgoooIACCrQiYIDVCqOZKKCAAgoooIACYwEDrLGFawoooIACCiigQCsCBlitMJqJAgoooIACCigwFjDAGlu4poACCiiggAIKtCJggNUKo5kooIACCiiggAJjAQOssYVrCiiggAIKKKBAKwIGWK0wmokCCiiggAIKKDAWMMAaW7imgAIKKKCAAgq0ImCA1QqjmSiggAIKKKCAAmMBA6yxhWsKKKCAAgoooEArAgZYrTCaiQIKKKCAAgooMBYwwBpbuKaAAgoooIACCrQiYIDVCqOZKKCAAgoooIACYwEDrLGFawoooIACCiigQCsCBlitMJqJAgoooIACCigwFjDAGlu4poACCiiggAIKtCJggNUKo5kooIACCiiggAJjAQOssYVrCiiggAIKKKBAKwIGWK0wmokCCiiggAIKKDAWMMAaW7imgAIKKKCAAgq0ImCA1QqjmSiggAIKKKCAAmMBA6yxhWsKKKCAAgoooEArAgZYrTCaiQIKKKCAAgooMBYwwBpbuKaAAgoooIACCrQiYIDVCqOZKKCAAgoooIACYwEDrLGFawoooIACCiigQCsCBlitMJqJAgoooIACCigwFjDAGlu4poACCiiggAIKtCJggNUKo5kooIACCiiggAJjAQOssYVrCiiggAIKKKBAKwIGWK0wmokCCiiggAIKKDAWMMAaW7imgAIKKKCAAgq0ImCA1QqjmSiggAIKKKCAAmMBA6yxhWsKKKCAAgoooEArAgZYrTCaiQIKKKCAAgooMBYwwBpbuKZA1wQGwDWBY7p2YB6PAgoo0LLAUcC128xzT5uZmZcCCiydwPWBRwDfAz4NfBb47tLtpTukgAIK7K7ANYBHAyl4+kz1+Aawf6u7ZYC1VTk/p8DyC5QqqLo58Dwgz78O/DPwAeCfqmBr33ZuIsvP4B4qoIACGwp8DTgT+FPgydW98gvAB6vXc+/MvfLyDXOqNpgMsHIDbqQEbnPn1fjc8q3mwA5cvt3q/R4d1HuBLQFs5hfVZcBfAxcDLwBOA+4BPAY4H/gk8A/AR4BvAWcDl25pr/rzodxOJu6VOfgrlkpgD6SK2KTApgQSFAy6c+3k7zT3wHlT7oW/VgVZdwXuBDwcuABIsPVe4B+r4Osc4MfrZdwMsPLHeMjajfccBKdfBAdesKreB3HZ4FAuGeTAPsahBfZMuTGuPWqf7ZzAhRT2cfEBV+dyLmVPuXh4jvy/sP4ZuPGJ8KdnrL/N5Ls/3gN//HF41fFw9l7gatXjesBDgB9VReIfBj4FnAV8aaMbyOS39OR5fqulvcZEunMKBC+E0ryvTmyzM0/3ceAl76Zc/m32d+MX8s6w+S3A+zmYK7j8ErgiQcWKp2sfDc85A47PD8xNpJd8Al50C/ji1YHDq8d1gNOrUqzPA7lX5gfq56rg6weTX9D4T1Zyw/0EcOO1G+27BC6/bFQtufadVXh2EPsGh3DJMCC/gkNKYa8B1lKduP0cwEUHHDgMsPaWSzikrOq1tnOsBx4Aext/u/N+8+UFXn8w/PYe+M56H7oI+GoVYOUGkmLzBF0XdqZIe72j3/C9ciTwJuA+aze9bB9ceukyXL8DymAvpRwwbE6ydi99psB6AoUySD1YYbDi/ytT0L9nAHsP2EIB0X5431546t5RDDVTLKVjaaf15eoHaqoTPwp8P/fKxk16GGDlZnqjmVn5hgIKrLhAqrFeDzyDDYKs+jhTXZhfZv9W/WJ7X6M6sadVicN75ZuBe9VILhVQoIsC7weeWpU9bXh8ubmeB6TqMD9Iz9z1ouwNd9kNFFCgRYHUbqUq66bzBlipTjyueqSU+7rAsVWUlp6JJgUUUKCjAqcAp84bYOXmmntjHmlCcPxkgNUo0eqol4elQK8Fvg08u2qnuS5EqgdSU/CVquTq41Vbg7Q9OHfdT/bjTe+V/TjPHmVvBdIs9U+Av9lIIPfKlF59syrdTxVhemp/sRlgZYPUI6bisq57zTI9jTbZQGyj/fF9BRTYYYECX78GPOMkeNOhM3q85e88DVsTVKUqMEM5ZD09C3O3MY0E0nA898qU5jVT7NJOzaSAAqsrUOCHh8Mf3gT+8uoz/qTTPCJ/7/nFmvtkehfmx2fGGLyysXvjV1jJ+mHVaAZ1gBWirDefry6be65APwX2w/2vDx9+Hlzw01CaMzjkZlA30vzQsBPRKKi6ZJPdm3skO7xXpmNy8wdq7p/NH6c98vBQFeiMQIF/fxS8/ulw3pPgisREdcqQDP9SdQDKsDZpoJWSqvyo6ml71JrGpQL9FUhJy2urm0BKqtNg/R3Ac4FHAicDDknW3+vDI1dAgZFAhmX4o6rmLj+Y0mg9pfl/DjwOuC1wqFgKKKBABI4GXll1G34b8BTg7sCJ1ZQQKimggAIKjH5kPrNqEpHBRH8PuG81dNXBWwH6/6BmItAP5ljvAAAAAElFTkSuQmCC"}}},{"metadata":{"trusted":true},"cell_type":"code","source":"normalized_patient_100 = histogram_normalize_4d(rescaled_patient_100)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_patient_slices_3d(normalized_patient_100[5])\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"### 3. Crop images aroind ROI \n\nIdentified using Fourier Transform over time"},{"metadata":{"trusted":true},"cell_type":"code","source":"# This is set after experiments. Could be lowered for shorter runtime\nheart_pixel_size = 150","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Fourier transform calculating first harmonic mean"},{"metadata":{"trusted":true},"cell_type":"code","source":"ft_patient_100 = fourier_time_transform(normalized_patient_100)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_patient_slices_3d(ft_patient_100)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Threshold for segmentation + erode for removing \"noise\" + mean coordinates for labels\n# Segmentation is not necessary here? We can erode than just calculate mean coordinates for intensity. \n# This would yeild a better estimate of ROI because the higher intensity the higher movement.\n\ny, x = roi_mean_yx(ft_patient_100)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"num_slices, time, h, w = normalized_patient_100.shape\nheart_cropped_patient_100_4d = np.zeros((num_slices, time, heart_pixel_size, heart_pixel_size))\n\nfor i in range(num_slices):\n        for j in range(time):\n            heart_cropped_patient_100_4d[i,j] = crop_roi(normalized_patient_100[i,j], heart_pixel_size, heart_pixel_size, y, x)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_patient_slices_3d(heart_cropped_patient_100_4d[5])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 4. Rotate images \nSuch that left ventricle is in top right part of img\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"rotated_patient_100 = rotate_images_210_deg(heart_cropped_patient_100_4d, patient_100.orientation)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_patient_slices_3d(rotated_patient_100[5])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 5. Segment out left ventricle\nFor each 2d slice. Uses watershed (go up for implementation). \nAlso tried Threshold, K-means and EM Segmentation with worse result.\nEM segmentation took way to long time. \n\nK-means similar result but way longer time."},{"metadata":{"trusted":true},"cell_type":"code","source":"segmented_left_ventricle_4d, labels = find_left_ventricle(rotated_patient_100,kmeans=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"kmeans_segmented_left_ventricle_4d, labels = find_left_ventricle(rotated_patient_100,kmeans=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Threshold Segmentation Plot"},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_patient_slices_3d(segmented_left_ventricle_4d[5])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Kmeans Segmentation Plot"},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_patient_slices_3d(kmeans_segmented_left_ventricle_4d[5])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 6. Caluculating Volume\nSince every 4d image is scaled 1:1 we sum upp all voxels for each timestep. Then we take max/min to calculate ejection rate."},{"metadata":{"trusted":true},"cell_type":"code","source":"def 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    return volume ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"v = volume_for_patient(segmented_left_ventricle_4d, patient_100.dist)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.plot(v)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"Just one bad segmentation ca lead to patient beeing miss diagnosed since max/min value where taken."},{"metadata":{},"cell_type":"markdown","source":"### Example with other patients"},{"metadata":{"trusted":true},"cell_type":"code","source":"patient_ex = load_patient(20)\nrescaled_patient_ex = rescale_patient_4d_imgs(patient_ex)\nnormalized_patient_ex = histogram_normalize_4d(rescaled_patient_ex)\nplot_patient_slices_3d(normalized_patient_ex[5])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ft_patient_ex = fourier_time_transform(normalized_patient_ex)\ny, x = roi_mean_yx(ft_patient_ex)\n\nnum_slices, time, h, w = normalized_patient_ex.shape\nheart_cropped_patient_ex_4d = np.zeros((num_slices, time, heart_pixel_size, heart_pixel_size))\n\nfor i in range(num_slices):\n        for j in range(time):\n            heart_cropped_patient_ex_4d[i,j] = crop_roi(normalized_patient_ex[i,j], heart_pixel_size, heart_pixel_size, y, x)\n            \nrotated_patient_ex = rotate_images_210_deg(heart_cropped_patient_ex_4d, patient_ex.orientation)\nsegmented_left_ventricle_4d_ex, labels = find_left_ventricle(rotated_patient_ex,kmeans=False)\nplot_patient_slices_3d(segmented_left_ventricle_4d_ex[5])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"patient_ex = load_patient(400)\nrescaled_patient_ex = rescale_patient_4d_imgs(patient_ex)\nnormalized_patient_ex = histogram_normalize_4d(rescaled_patient_ex)\nplot_patient_slices_3d(normalized_patient_ex[5])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ft_patient_ex = fourier_time_transform(normalized_patient_ex)\ny, x = roi_mean_yx(ft_patient_ex)\n\nnum_slices, time, h, w = normalized_patient_ex.shape\nheart_cropped_patient_ex_4d = np.zeros((num_slices, time, heart_pixel_size, heart_pixel_size))\n\nfor i in range(num_slices):\n        for j in range(time):\n            heart_cropped_patient_ex_4d[i,j] = crop_roi(normalized_patient_ex[i,j], heart_pixel_size, heart_pixel_size, y, x)\n            \nrotated_patient_ex = rotate_images_210_deg(heart_cropped_patient_ex_4d, patient_ex.orientation)\nsegmented_left_ventricle_4d_ex, labels = find_left_ventricle(rotated_patient_ex,kmeans=False)\nplot_patient_slices_3d(segmented_left_ventricle_4d_ex[5])","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}