From 1920c99bd2ce8cd124d64a5cbb6c44ceea9967ae Mon Sep 17 00:00:00 2001 From: s224362 <s224362@student.dtu.dk> Date: Tue, 7 Jan 2025 12:57:39 +0100 Subject: [PATCH] Added some notebooks to try using gitlab --- .../Agamodon_anguliceps_processing.ipynb | 1089 +++++++++++ Notebooks/Beanies_processing.ipynb | 1726 +++++++++++++++++ Notebooks/Curve_Analysis_Script.ipynb | 380 ++++ 3 files changed, 3195 insertions(+) create mode 100644 Notebooks/Agamodon_anguliceps_processing.ipynb create mode 100644 Notebooks/Beanies_processing.ipynb create mode 100644 Notebooks/Curve_Analysis_Script.ipynb diff --git a/Notebooks/Agamodon_anguliceps_processing.ipynb b/Notebooks/Agamodon_anguliceps_processing.ipynb new file mode 100644 index 0000000..0ec2772 --- /dev/null +++ b/Notebooks/Agamodon_anguliceps_processing.ipynb @@ -0,0 +1,1089 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Imports" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from mpl_toolkits.mplot3d.art3d import Poly3DCollection\n", + "from skimage import measure\n", + "from scipy.interpolate import RegularGridInterpolator\n", + "from scipy import signal\n", + "import open3d as o3d\n", + "import sympy as sp\n", + "from scipy.optimize import curve_fit\n", + "import sympy as sp\n", + "from mpl_toolkits.mplot3d import Axes3D\n", + "import nibabel as nib\n", + "import matplotlib.animation as animation\n", + "from IPython.display import HTML\n", + "import matplotlib\n", + "from ipywidgets import interactive, fixed, IntSlider\n", + "import math\n", + "from skimage.morphology import dilation, ball\n", + "from skimage.morphology import skeletonize_3d, skeletonize\n", + "from matplotlib.widgets import RectangleSelector\n", + "from scipy.ndimage import gaussian_filter\n", + "from mayavi import mlab\n", + "import nibabel as nib\n", + "from matplotlib.ticker import MaxNLocator\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Helper Functions" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Function to marching cube meshes\n", + "def plotMesh(ax, mesh, ax_x, ax_y, ax_z, azim, elev):\n", + " \n", + " ax.add_collection3d(mesh)\n", + " ax.set_xlabel(\"x\")\n", + " ax.set_ylabel(\"y\")\n", + " ax.set_zlabel(\"z\")\n", + " ax.set_xlim(ax_x[0], ax_x[1]) \n", + " ax.set_ylim(ax_y[0], ax_y[1]) \n", + " ax.set_zlim(ax_z[0], ax_z[1]) \n", + " ax.set_aspect('equal')\n", + " ax.azim = azim\n", + " ax.elev = elev\n", + "\n", + "# Define the 3D polynomial function of 4th degree\n", + "def poly_3d(xy, a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13, a14, a15):\n", + " x, y = xy\n", + " return (a0 + a1*x + a2*y + a3*x**2 + a4*y**2 + a5*x*y + a6*x**3 + a7*y**3 + a8*x**2*y + a9*x*y**2 +\n", + " a10*x**4 + a11*y**4 + a12*x**3*y + a13*x**2*y**2 + a14*x*y**3)\n", + "\n", + "# Fits polynomial to set of datapoints and displays normalvectors (plot)\n", + "def fit_polynomial_surface(points, subset=10, normal_direction=0, num_points_quiver=10, num_points_surface=50, scale_factor=35):\n", + " '''\n", + " num_points_quiver = 10 # Adjust the number of points as needed for the quivers\n", + " num_points_surface = 50 # Adjust the number of points as needed for the surface\n", + " scale_factor = 35 # Adjust the scale of the normal vectors\n", + " normal_direction = 1 # Change this to 1 (or 0) to reverse the direction of the normals\n", + " WANT NORMALS FACING INTO BEANIE SKULL\n", + " \n", + " extra_distance: The extra distance added around the points for sampling\n", + " step_size: The step size between each sampling along the normal vectors\n", + " steps: The amount of steps done for sampling along normal vectors\n", + " grid_size_x: Size of the grid of sample points along the x-axis (can be left out and calculated automatically)\n", + " grid_size_y: Size of the grid of sample points along the y-axis (can be left out and calculated automatically)\n", + " '''\n", + "\n", + " # Extract x, y, z coordinates from the points (These are subsets defined earlier!)\n", + " x_num = points[::subset, 0]\n", + " y_num = points[::subset, 1]\n", + " z_num = points[::subset, 2]\n", + "\n", + " # Initial guess for the parameters\n", + " p0 = np.ones(16) # Needs to match the amount of a's in poly_3d (helper function)\n", + "\n", + " # Use curve_fit to fit the function to the data\n", + " popt, pcov = curve_fit(poly_3d, (x_num, y_num), z_num, p0)\n", + "\n", + " # Create a grid of x, y values for the surface\n", + " x_grid_surface, y_grid_surface = np.meshgrid(np.linspace(min(x_num), max(x_num), num_points_surface), np.linspace(min(y_num), max(y_num), num_points_surface))\n", + "\n", + " # Compute the corresponding z values for the surface\n", + " z_grid_surface = poly_3d((x_grid_surface, y_grid_surface), *popt)\n", + "\n", + " # Create a grid of x, y values for the quivers\n", + " x_grid_quiver, y_grid_quiver = np.meshgrid(np.linspace(min(x_num), max(x_num), num_points_quiver), np.linspace(min(y_num), max(y_num), num_points_quiver))\n", + "\n", + " # Compute the corresponding z values for the quivers\n", + " z_grid_quiver = poly_3d((x_grid_quiver, y_grid_quiver), *popt)\n", + "\n", + " # Define the symbols\n", + " a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13, a14, a15 = sp.symbols('a0 a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 a11 a12 a13 a14 a15')\n", + " x, y = sp.symbols('x y')\n", + "\n", + " # Define the polynomial function\n", + " poly_expr = (a0 + a1*x + a2*y + a3*x**2 + a4*y**2 + a5*x*y + a6*x**3 + a7*y**3 + a8*x**2*y + a9*x*y**2 +\n", + " a10*x**4 + a11*y**4 + a12*x**3*y + a13*x**2*y**2 + a14*x*y**3)\n", + "\n", + " # Differentiate with respect to x and y\n", + " dz_dx = sp.diff(poly_expr, x)\n", + " dz_dy = sp.diff(poly_expr, y)\n", + "\n", + " # Create lambdified functions for evaluating the derivatives\n", + " dz_dx_func = sp.lambdify((x, y, a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13, a14, a15), dz_dx)\n", + " dz_dy_func = sp.lambdify((x, y, a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13, a14, a15), dz_dy)\n", + "\n", + " # Evaluate the derivatives at each point on the quiver surface\n", + " dz_dx_vals = dz_dx_func(x_grid_quiver, y_grid_quiver, *popt)\n", + " dz_dy_vals = dz_dy_func(x_grid_quiver, y_grid_quiver, *popt)\n", + "\n", + " # Compute the normal vectors\n", + " normals = np.dstack((-dz_dx_vals, -dz_dy_vals, np.ones_like(dz_dx_vals)))\n", + "\n", + " # To change direction of normals\n", + " normals *= (-1)**normal_direction \n", + "\n", + " # Normalize the normals\n", + " normals /= np.linalg.norm(normals, axis=2, keepdims=True)\n", + "\n", + " # Create a 3D plot\n", + " fig = plt.figure()\n", + " ax = fig.add_subplot(111, projection='3d')\n", + "\n", + " # Plot the original points\n", + " ax.scatter(x_num, y_num, z_num, color='b')\n", + "\n", + " # Plot the fitted surface using the surface grid\n", + " ax.plot_surface(x_grid_surface, y_grid_surface, z_grid_surface, color='r', alpha=0.5)\n", + "\n", + " # Plot the normal vectors using the quiver grid\n", + " for i in range(len(x_grid_quiver)):\n", + " for j in range(len(y_grid_quiver)):\n", + " ax.quiver(x_grid_quiver[i, j], y_grid_quiver[i, j], z_grid_quiver[i, j], normals[i, j, 0], normals[i, j, 1], normals[i, j, 2],\n", + " length=scale_factor, linewidth=5, color='g')\n", + "\n", + " ax.set_aspect('equal')\n", + " plt.show()\n", + " \n", + " return popt, dz_dx_func, dz_dy_func\n", + "\n", + "# Generates PCA-coordinates to sample \n", + "def generate_sample_points(points, popt, dz_dx_func, dz_dy_func, normal_direction=0, subset=10, extra_distance=15, step_size=0.1, min_steps=10, grid_size_x=None, grid_size_y=None):\n", + " '''\n", + " subset: amount of subsamples used\n", + " extra_distance: The extra distance added around the points for sampling\n", + " step_size: The step size between each sampling along the normal vectors\n", + " steps: The amount of steps done for sampling along normal vectors\n", + " grid_size_x: Size of the grid of sample points along the x-axis (can be left out and calculated automatically)\n", + " grid_size_y: Size of the grid of sample points along the y-axis (can be left out and calculated automatically)\n", + "\n", + " STANDARD FOR AGAMODON:\n", + " extra_distance = 15\n", + " step_size = 0.1\n", + " min_steps: The number of steps in the least sampled direction! (The other direction is sampled 5x as much)\n", + " grid_size_x = 500\n", + " grid_size_y = 200\n", + " '''\n", + " # Extract x, y, z coordinates from the points (These are subsets defined earlier!)\n", + " x_num = points[::subset, 0]\n", + " y_num = points[::subset, 1]\n", + " z_num = points[::subset, 2]\n", + "\n", + " # For calculating the span of the points to generate meshgrid\n", + " min_values = np.min(points, axis=0)\n", + " max_values = np.max(points, axis=0)\n", + "\n", + " # If user has not chosen, make it the size of the spans\n", + " if grid_size_x is None:\n", + " grid_size_x = int(max_values[0] - min_values[0])\n", + " if grid_size_y is None:\n", + " grid_size_y = int(max_values[1] - min_values[1])\n", + " \n", + " # Step 1: Create a grid in the x-y-plane\n", + " x_grid, y_grid = np.meshgrid(np.linspace(min(x_num) - extra_distance, max(x_num) + extra_distance, grid_size_x), np.linspace(min(y_num) - extra_distance, max(y_num) + extra_distance, grid_size_y))\n", + "\n", + " # Step 2: Compute the height of the fitted polynomial surface at each point\n", + " z_grid = poly_3d((x_grid, y_grid), *popt)\n", + "\n", + " # Step 3: Compute the normal vector at each point on the grid\n", + " dz_dx_vals = dz_dx_func(x_grid, y_grid, *popt) # defined priorly. \n", + " dz_dy_vals = dz_dy_func(x_grid, y_grid, *popt)\n", + " normals = np.dstack((-dz_dx_vals, -dz_dy_vals, np.ones_like(dz_dx_vals)))\n", + "\n", + " normals *= (-1)**normal_direction # To change direction of normals\n", + "\n", + " # Step 4: Normalize the normals for better visualization\n", + " normals /= np.linalg.norm(normals, axis=2, keepdims=True)\n", + "\n", + " # Step 5: Go `steps` steps in the direction of the normal vector and `steps` steps in the opposite direction to create the layers\n", + " # Modify the number of steps in each direction\n", + " steps_positive = 5 * min_steps\n", + " steps_negative = min_steps\n", + " grid_3d = np.zeros((grid_size_y, grid_size_x, steps_positive + steps_negative + 1, 3)) # 3D array to store the x, y, and z coordinates of the points (each index in the 3D matrix contains a 3D-coordinate in PCA-space)\n", + "\n", + " # Go steps_negative steps in the opposite direction\n", + " for i in range(-steps_negative, 0):\n", + " displacement = i * step_size * normals # displacement in x, y, and z directions\n", + " grid_3d[:, :, i+steps_negative, 0] = x_grid + displacement[:, :, 0] # x coordinate\n", + " grid_3d[:, :, i+steps_negative, 1] = y_grid + displacement[:, :, 1] # y coordinate\n", + " grid_3d[:, :, i+steps_negative, 2] = z_grid + displacement[:, :, 2] # z coordinate\n", + "\n", + " # Go steps_positive steps in the direction of the normal vector\n", + " for i in range(steps_positive):\n", + " displacement = i * step_size * normals # displacement in x, y, and z directions\n", + " grid_3d[:, :, i+steps_negative+1, 0] = x_grid + displacement[:, :, 0] # x coordinate\n", + " grid_3d[:, :, i+steps_negative+1, 1] = y_grid + displacement[:, :, 1] # y coordinate\n", + " grid_3d[:, :, i+steps_negative+1, 2] = z_grid + displacement[:, :, 2] # z coordinate\n", + "\n", + " # Flatten the 3D grid\n", + " points_flat = grid_3d.reshape(-1, 3)\n", + "\n", + " # Plot a small subset of the found points\n", + " fig = plt.figure()\n", + " ax = fig.add_subplot(111, projection='3d')\n", + " ax.scatter(points_flat[::2000, 0], points_flat[::2000, 1], points_flat[::2000, 2])\n", + " plt.show()\n", + "\n", + " return points_flat, grid_3d.shape[0:3]\n", + "\n", + "# Interpolate sample points in original volume file\n", + "def sample_in_original_volume(points_original, original_nii_path, grid_shape, intMethod='linear', expVal=0.0):\n", + " \n", + " # Load in original volume\n", + " original = nib.load(original_nii_path)\n", + " imgDim = original.header['dim'][1:4]\n", + "\n", + " # Set-up interpolators for moving image\n", + " x = np.arange(start=0, stop=imgDim[0], step=1)\n", + " y = np.arange(start=0, stop=imgDim[1], step=1)\n", + " z = np.arange(start=0, stop=imgDim[2], step=1)\n", + " F_moving = RegularGridInterpolator((x, y, z), original.get_fdata().astype('float16'), method=intMethod, bounds_error=False, fill_value=expVal)\n", + "\n", + " # Evaluate transformed grid points in the moving image\n", + " fVal = F_moving(points_original)\n", + " # Reshape to voxel grid\n", + " volQ = np.reshape(fVal,newshape=grid_shape).astype('float32')\n", + "\n", + " return volQ\n", + "\n", + "# Animate moving through the layers of volume\n", + "def ani_through_volume(volume, indices=None):\n", + " '''\n", + " indices: List [a,b] of layer numbers you want animated\n", + " '''\n", + " def update_img(z):\n", + " ax.clear()\n", + " ax.imshow(volume[:,:,z])\n", + " ax.set_title('Layer ' + str(z))\n", + "\n", + " # User can specify indices they want to see\n", + " if indices != None:\n", + " a = indices[0]\n", + " b = indices[1]\n", + " else:\n", + " a = 0\n", + " b = volume.shape[2] - 1\n", + "\n", + " matplotlib.rcParams['animation.embed_limit'] = 2**128\n", + "\n", + " fig, ax = plt.subplots()\n", + " ani = animation.FuncAnimation(fig, update_img, frames=range(a,b+1), interval=200)\n", + " \n", + " return HTML(ani.to_jshtml())\n", + "\n", + "# Create interactive plot for slicing\n", + "def interactive_plot_slicing(volume):\n", + " def get_slice(vol, x_start, x_end, y_start, y_end, z):\n", + " slice = vol[x_start:x_end, y_start:y_end, z]\n", + " plt.imshow(slice)\n", + " plt.show()\n", + "\n", + " x_start_slider = IntSlider(min=0, max=volume.shape[0]-1, value=0)\n", + " x_end_slider = IntSlider(min=0, max=volume.shape[0]-1, value=volume.shape[0]-1)\n", + " y_start_slider = IntSlider(min=0, max=volume.shape[1]-1, value=0)\n", + " y_end_slider = IntSlider(min=0, max=volume.shape[1]-1, value=volume.shape[1]-1)\n", + " z_slider = IntSlider(min=0, max=volume.shape[2]-1, value=0)\n", + "\n", + " def update_x_end_range(*args):\n", + " x_end_slider.min = x_start_slider.value\n", + " x_start_slider.observe(update_x_end_range, 'value')\n", + "\n", + " def update_y_end_range(*args):\n", + " y_end_slider.min = y_start_slider.value\n", + " y_start_slider.observe(update_y_end_range, 'value')\n", + "\n", + " interactive_plot = interactive(get_slice, \n", + " vol=fixed(volume), \n", + " x_start=x_start_slider, \n", + " x_end=x_end_slider, \n", + " y_start=y_start_slider, \n", + " y_end=y_end_slider, \n", + " z=z_slider)\n", + " return interactive_plot\n", + "\n", + "# Allows user to remove unwanted sutures\n", + "def interactive_process_volume(volume, interactive_plot, z_vals):\n", + " # Define the slicing parameters\n", + " x_start = interactive_plot.kwargs['x_start']\n", + " x_end = interactive_plot.kwargs['x_end']\n", + " y_start = interactive_plot.kwargs['y_start']\n", + " y_end = interactive_plot.kwargs['y_end']\n", + " \n", + " min_z, max_z = z_vals[0], z_vals[1]\n", + "\n", + " # Crop the volume based on the provided slices\n", + " cropped_volume = volume[x_start:x_end, y_start:y_end, min_z:max_z+1]\n", + "\n", + " # Calculate the median layer index\n", + " median_z = cropped_volume.shape[2] // 2\n", + " median_layer = cropped_volume[:, :, median_z]\n", + "\n", + " # Create a new volume to apply changes\n", + " new_volume = cropped_volume.copy()\n", + "\n", + " fig, ax = plt.subplots()\n", + " ax.imshow(median_layer, cmap='gray')\n", + "\n", + " def line_select_callback(eclick, erelease):\n", + " nonlocal new_volume\n", + " x1, y1 = int(eclick.xdata), int(eclick.ydata)\n", + " x2, y2 = int(erelease.xdata), int(erelease.ydata)\n", + " if x1 > x2:\n", + " x1, x2 = x2, x1\n", + " if y1 > y2:\n", + " y1, y2 = y2, y1\n", + "\n", + " # Find the maximum value within the selected area\n", + " max_value = np.max(new_volume[y1:y2+1, x1:x2+1, :])\n", + "\n", + " # Assign the maximum value to all pixels within the selected area for each slice\n", + " for z in range(new_volume.shape[2]):\n", + " new_volume[y1:y2+1, x1:x2+1, z] = max_value\n", + "\n", + " # Update the displayed image in real-time\n", + " ax.clear()\n", + " ax.imshow(new_volume[:, :, median_z], cmap='gray')\n", + " plt.draw()\n", + "\n", + " rs = RectangleSelector(ax, line_select_callback,\n", + " useblit=True,\n", + " button=[1], # Only use left button\n", + " minspanx=5, minspany=5,\n", + " spancoords='pixels',\n", + " interactive=True)\n", + "\n", + " plt.show(block=True) # Block until plot window is closed\n", + "\n", + " # Return the new_volume after the plot is closed\n", + " return new_volume\n", + "\n", + "# Interactive plot for deciding threshold value\n", + "def interactive_plot_threshold(volume):\n", + "\n", + " # Define the function to generate the plot\n", + " def plot_threshold(threshold):\n", + " slice = volume[:,:, volume.shape[2]//2]\n", + " slice_bin = np.where(slice >= threshold, 1, 0)\n", + " slice_bin= 1 - slice_bin # Make suture the \"foreground\"\n", + " plt.imshow(slice_bin)\n", + " plt.show()\n", + "\n", + " # Calculate the minimum value for the slider\n", + " min_val = (volume.max()+volume.min())//2\n", + " # Round up to the nearest hundred\n", + " min_val = math.ceil(min_val / 100.0) * 100\n", + "\n", + " # Create a slider for the threshold value\n", + " threshold_slider = IntSlider(min=min_val, max=volume.max(), step=100, value=min_val)\n", + "\n", + " # Create the interactive plot\n", + " interactive_plot_new = interactive(plot_threshold, threshold=threshold_slider)\n", + "\n", + " return interactive_plot_new\n", + "\n", + "# Extracts the suture depending on your former choices\n", + "def extract_suture(volume, points_orig, interactive_plot_slice, interactive_plot_thresh, z_vals, grid_shape):\n", + " '''\n", + " Input volume is sliced\n", + " '''\n", + " min_z = z_vals[0]\n", + " max_z = z_vals[1]\n", + " num_indices = max_z - min_z\n", + "\n", + " # Extrac chosen indices\n", + " x_start = interactive_plot_slice.kwargs['x_start']\n", + " x_end = interactive_plot_slice.kwargs['x_end']\n", + " y_start = interactive_plot_slice.kwargs['y_start']\n", + " y_end = interactive_plot_slice.kwargs['y_end']\n", + "\n", + " # Extract chosen threshold\n", + " threshold_value = interactive_plot_thresh.kwargs['threshold']\n", + " \n", + " # Generate a slice\n", + " slice = volume[:,:, volume.shape[2]//2]\n", + "\n", + " # Generate empty array to store suture points in\n", + " suture = np.zeros((slice.shape[0], slice.shape[1], volume.shape[2]))\n", + " # For animation\n", + " binary_slices = np.zeros((slice.shape[0], slice.shape[1], num_indices))\n", + "\n", + " # Two first indicies were bad in segmentation\n", + " for index in range(num_indices):\n", + " slice = volume[:,:, index]\n", + " slice_bin = np.where(slice >= threshold_value, 1, 0)\n", + " slice_bin= 1 - slice_bin # Make suture the \"foreground\"\n", + "\n", + " # Label the connected components in the binary image\n", + " labels = measure.label(slice_bin)\n", + "\n", + " # Now, you can analyze each individual blob\n", + " properties = measure.regionprops(labels)\n", + "\n", + " # Sort the regions by area_bbox in descending order\n", + " properties_sorted = sorted(properties, key=lambda x: x.bbox_area, reverse=True)\n", + "\n", + " # Get the region with the largest area_bbox\n", + " largest_bbox_region = properties_sorted[0]\n", + "\n", + " # Get the coordinates of the pixels in the blob with the largest area_bbox\n", + " coords_largest_bbox = largest_bbox_region.coords\n", + " \n", + " # Extract x and y coordinates\n", + " x_coords = coords_largest_bbox[:, 1]\n", + " y_coords = coords_largest_bbox[:, 0]\n", + "\n", + " suture[y_coords, x_coords, index] = 1\n", + " binary_slices[:,:, index] = slice_bin\n", + "\n", + " ### ANIMATION\n", + " # Create a figure with two subplots\n", + " fig, axs = plt.subplots(1, 2, figsize=(8,3))\n", + "\n", + " # Initialize the images\n", + " im1 = axs[0].imshow(suture[:, :, 0], animated=True)\n", + " im2 = axs[1].imshow(volume[:,:, 0], animated=True)\n", + "\n", + " # Initialize the text\n", + " txt = axs[0].text(0.5, 1.01, f'Layer: {0}', transform=axs[0].transAxes)\n", + "\n", + " # Update function for the animation\n", + " def updatefig(i):\n", + " im1.set_array(suture[:, :, i])\n", + " im2.set_array(volume[:,:, i])\n", + " txt.set_text(f'Layer: {i}')\n", + " return im1, im2, txt\n", + "\n", + " # Create the animation\n", + " ani = animation.FuncAnimation(fig, updatefig, frames=range(num_indices), interval=200, blit=True)\n", + " \n", + " ### GET ORIGINAL COORDINATES\n", + " new_shape = grid_shape + (3,)\n", + " coords_orig = points_orig.reshape(new_shape)\n", + "\n", + " # cooresponding original coordinates to sliced area\n", + " coords_orig_sliced = coords_orig[x_start:x_end, y_start:y_end, min_z:min_z+num_indices]\n", + " \n", + " ### PLOT ORIGINAL COORDINATES OF SLICED AREA\n", + " # Get the indices where suture is 1\n", + " indices = np.where(suture == 1)\n", + "\n", + " # Get the corresponding points in coords_orig_sliced\n", + " points = coords_orig_sliced[indices]\n", + "\n", + " # Create a 3D scatter plot\n", + " fig = plt.figure()\n", + " ax = fig.add_subplot(111, projection='3d')\n", + " ax.scatter(points[::10, 0], points[::10, 1], points[::10, 2])\n", + " ax.set_xlabel('x')\n", + " ax.set_ylabel('y')\n", + " ax.set_zlabel('z')\n", + " plt.show()\n", + "\n", + " return ani, points, suture, coords_orig_sliced\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## \"unfolded\" volume generation" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plot point-cloud from input file" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "input_file = 'Agamodon_anguliceps_cut_PLY.ply'\n", + "\n", + "pcd = o3d.io.read_point_cloud(input_file) \n", + "points = np.asarray(pcd.points)\n", + "\n", + "\n", + "# Plot a subset of the points\n", + "# Choose a subset of points and view point of display\n", + "subset = 10\n", + "viewPoint = [90, -45]\n", + "\n", + "fig = plt.figure(figsize=(6, 5))\n", + "ax = fig.add_subplot(1,1,1, projection='3d')\n", + "ax.scatter(points[::subset,0], points[::subset,1], points[::subset,2], c='k', marker='.')\n", + "ax.set_xlabel(\"x\")\n", + "ax.set_ylabel(\"y\")\n", + "ax.set_zlabel(\"z\")\n", + "ax.set_aspect('auto')\n", + "ax.azim = viewPoint[0]\n", + "ax.elev = viewPoint[1]\n", + "ax.set_title('Mesh-points in original domain')\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Transform point-cloud into PCA-space" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Principal Component Analysis (PCA) of the point cloud\n", + "points_mu = np.mean(points, axis = 0) # Center the point cloud\n", + "cov = np.cov(np.transpose(points)) # Covariance matrix \n", + "\n", + "# SVD or Eigendecomposition\n", + "U,S,V = np.linalg.svd(cov,full_matrices=False)\n", + "\n", + "# Apply transform to rotated pointcloud\n", + "points_rot = (points - points_mu) @ U" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Display points in PCA-space" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib qt5" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Plot the rotated pointcloud\n", + "fig = plt.figure(figsize=(8,4))\n", + "\n", + "# display the original point cloud with the principal axes\n", + "ax = fig.add_subplot(1,2,1, projection='3d')\n", + "ax.scatter(points[::subset,0], points[::subset,1], points[::subset,2], c='k', marker='.')\n", + "# Set the number of ticks on each axis\n", + "ax.xaxis.set_major_locator(MaxNLocator(nbins=5))\n", + "ax.yaxis.set_major_locator(MaxNLocator(nbins=5))\n", + "ax.zaxis.set_major_locator(MaxNLocator(nbins=5))\n", + "ax.set_xlabel(\"x\")\n", + "ax.set_ylabel(\"y\")\n", + "ax.set_zlabel(\"z\")\n", + "ax.set_aspect('equal')\n", + "ax.azim = 40\n", + "ax.elev = -25\n", + "ax.set_title('Original point cloud')\n", + "\n", + "# Add the principal axes to the plot using quiver\n", + "ax.quiver(points_mu[0], points_mu[1], points_mu[2], U[0,0], U[1,0], U[2,0], color='r', length=100, normalize=True, linewidth=2.5)\n", + "ax.quiver(points_mu[0], points_mu[1], points_mu[2], U[0,1], U[1,1], U[2,1], color='g', length=100, normalize=True, linewidth=2.5)\n", + "ax.quiver(points_mu[0], points_mu[1], points_mu[2], U[0,2], U[1,2], U[2,2], color='b', length=100, normalize=True, linewidth=2.5)\n", + "\n", + "# Display the rotated point cloud \n", + "ax = fig.add_subplot(1,2,2, projection='3d')\n", + "ax.scatter(points_rot[::subset,0], points_rot[::subset,1], points_rot[::subset,2], c='k', marker='.')\n", + "# Set the number of ticks on each axis\n", + "ax.xaxis.set_major_locator(MaxNLocator(nbins=5))\n", + "ax.yaxis.set_major_locator(MaxNLocator(nbins=5))\n", + "ax.zaxis.set_major_locator(MaxNLocator(nbins=2))\n", + "ax.set_xlabel(\"PC1\")\n", + "ax.set_ylabel(\"PC2\")\n", + "ax.set_zlabel(\"PC3\")\n", + "ax.set_aspect('equal')\n", + "ax.azim = 90\n", + "ax.elev = -45\n", + "ax.set_title('Rotated point cloud')\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Fit polynomial surface in PCA-space" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "\n", + "normal_direction = 1\n", + "\n", + "popt, dz_dx_func, dz_dy_func = fit_polynomial_surface(points_rot, normal_direction=normal_direction)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Generate sample point coordinates in PCA-space (and plot the sample points in PCA-space)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "points_flat, grid_shape = generate_sample_points(points_rot, popt, dz_dx_func, dz_dy_func, normal_direction=normal_direction, extra_distance=15, min_steps=20, step_size=0.4)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Transform sample points into world coordinates" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# From SVD earlier\n", + "points_orig = points_flat @ U.T + points_mu" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Visualise the sample points and the original surface" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib qt5" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "ax = fig.add_subplot(111, projection='3d')\n", + "ax.scatter(points_orig[::1000, 0], points_orig[::1000, 1], points_orig[::1000, 2], alpha=0.3)\n", + "ax.scatter(points[::30,0], points[::30,1], points[::30,2])\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Sample points_orig's values in the original volume" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "path_nii = 'G:/Mit drev/DTU/Fagprojekt/DATA/Agamodon_anguliceps_CAS153467/Agamodon_anguliceps_CAS153467000.nii'\n", + "volQ = sample_in_original_volume(points_orig, path_nii, grid_shape) # Can take some seconds" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Save the new sampled volume\n", + "niiPCA = nib.Nifti1Image(volQ, np.eye(4))\n", + "name = 'Agamodon_anguliceps_unfolded_volume_new.nii'\n", + "nib.save(niiPCA, name)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Suture extraction" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Load the NIfTI-file back in" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Load the .nii file in\n", + "name = 'Agamodon_anguliceps_unfolded_volume_new.nii'\n", + "nii_img = nib.load(name)\n", + "# Get the data as a numpy array\n", + "volQ = nii_img.get_fdata()\n", + "\n", + "# Print max and min values\n", + "print('Volume minimum: ', volQ.min())\n", + "print('Volume maximum: ', volQ.max())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Animate moving through the layers of the \"unfolded\" volume" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ani_through_volume(volQ)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Crop slices" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "interactive_plot_slice = interactive_plot_slicing(volQ)\n", + "interactive_plot_slice" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Remove unwanted sutures interactively" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib qt5\n", + "\n", + "z_vals = [60,61] # Choose z_vals you want to work on\n", + "\n", + "new_volume = interactive_process_volume(volQ, interactive_plot_slice, z_vals)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Choose layers and thresholding value" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "interactive_plot_thresh = interactive_plot_threshold(new_volume)\n", + "interactive_plot_thresh" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib qt5\n", + "ani, points, suture, coords_orig = extract_suture(new_volume, points_orig, interactive_plot_slice, interactive_plot_thresh, z_vals, grid_shape)\n", + "HTML(ani.to_jshtml())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Skeletonize" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "- Chooses only ONE layer to skeletonize and represent the suture (curve)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Choose layer\n", + "z_val = 0\n", + "\n", + "# Create skeleton\n", + "skeleton = skeletonize(suture[:,:,z_val])\n", + "\n", + "# Get the 2D coordinates of the skeleton points\n", + "y, x = np.where(skeleton)\n", + "\n", + "# Create a 2D plot\n", + "fig, ax = plt.subplots()\n", + "ax.scatter(x, y, color='r')\n", + "\n", + "# Set labels and title\n", + "ax.set_xlabel('X')\n", + "ax.set_ylabel('Y')\n", + "ax.set_title('2D Skeleton')\n", + "\n", + "# Show the plot\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Get the 3D coordinates of the skeleton points\n", + "skeleton_indices = np.where(skeleton == 1)\n", + "skeleton_coords = coords_orig[skeleton_indices + (z_val,)]\n", + "skeleton_coords = np.squeeze(skeleton_coords)\n", + "\n", + "# Separate the x, y, and z coordinates\n", + "x, y, z = skeleton_coords[:, 0], skeleton_coords[:, 1], skeleton_coords[:, 2]\n", + "\n", + "# Create a 3D plot\n", + "fig = plt.figure()\n", + "ax = fig.add_subplot(111, projection='3d')\n", + "ax.scatter(x, y, z, color='r')\n", + "\n", + "# Set labels and title\n", + "ax.set_xlabel('X')\n", + "ax.set_ylabel('Y')\n", + "ax.set_zlabel('Z')\n", + "ax.set_title('3D Skeleton')\n", + "\n", + "# Show the plot\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Export the skeleton" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "np.save(\"Agamodon_anguliceps_skeleton.npy\", skeleton_coords)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plot of original cranium and extracted suture" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Plot the original point cloud with the extracted skeleton line\n", + "\n", + "input_file = 'Agamodon_anguliceps_cut_PLY.ply'\n", + "\n", + "pcd = o3d.io.read_point_cloud(input_file) \n", + "points_from_input = np.asarray(pcd.points)\n", + "\n", + "\n", + "fig = plt.figure(figsize=(8, 8))\n", + "ax = fig.add_subplot(111, projection='3d')\n", + "ax.scatter(points_from_input[::subset,0], points_from_input[::subset,1], points_from_input[::subset,2], c='k', marker='.', alpha=0.7) # Make points semi-transparent\n", + "ax.scatter(skeleton_coords[:, 0], skeleton_coords[:, 1], skeleton_coords[:, 2], color='r', s=50) # Add mean_points to the plot with larger size\n", + "ax.set_xlabel(\"x\")\n", + "ax.set_ylabel(\"y\")\n", + "ax.set_zlabel(\"z\")\n", + "ax.set_aspect('equal')\n", + "ax.azim = 40\n", + "ax.elev = -25\n", + "ax.set_title('Original point cloud')\n", + "\n", + "# Add the principal axes to the plot using quiver\n", + "ax.quiver(points_mu[0], points_mu[1], points_mu[2], U[0,0], U[1,0], U[2,0], color='r', length=2, normalize=True, linewidth=2.5)\n", + "ax.quiver(points_mu[0], points_mu[1], points_mu[2], U[0,1], U[1,1], U[2,1], color='g', length=2, normalize=True, linewidth=2.5)\n", + "ax.quiver(points_mu[0], points_mu[1], points_mu[2], U[0,2], U[1,2], U[2,2], color='b', length=2, normalize=True, linewidth=2.5)\n", + "\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Mayavi show cranium with extracted points" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Load the .nii file\n", + "img = nib.load('Agamodon_anguliceps_CAS153467000.nii')\n", + "data = img.get_fdata()\n", + "\n", + "# Display the 3D volume\n", + "mlab.contour3d(data)\n", + "\n", + "# WITH EXTRACTED POINTS\n", + "mlab.points3d(points[:, 0], points[:, 1], points[:, 2], color=(1, 0, 0), scale_factor=30)\n", + "\n", + "mlab.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Mayavi show cranium with skeleton" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Load the .nii file\n", + "img = nib.load('Agamodon_anguliceps_CAS153467000.nii')\n", + "data = img.get_fdata()\n", + "\n", + "# Display the 3D volume\n", + "mlab.contour3d(data)\n", + "\n", + "# WITH SKELETON-POINTS\n", + "mlab.points3d(skeleton_coords[:, 0], skeleton_coords[:, 1], skeleton_coords[:, 2], color=(1, 0, 0), scale_factor=32)\n", + "\n", + "mlab.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Display extracted points and skeleton in PCA-space" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from sklearn.decomposition import PCA\n", + "import matplotlib.pyplot as plt\n", + "\n", + "# Assuming your data is in a variable called data\n", + "pca = PCA(n_components=2)\n", + "pca_points = pca.fit_transform(points)\n", + "\n", + "# Plot the transformed points\n", + "plt.scatter(pca_points[:, 0], pca_points[:, 1])\n", + "\n", + "# Transform the skeleton_coords into the same PCA space\n", + "pca_skeleton_coords = pca.transform(skeleton_coords)\n", + "\n", + "# Plot the skeleton_coords in the same PCA space\n", + "plt.scatter(pca_skeleton_coords[:, 0], pca_skeleton_coords[:, 1], color='r')\n", + "\n", + "plt.show()" + ] + } + ], + "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.11.0" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/Notebooks/Beanies_processing.ipynb b/Notebooks/Beanies_processing.ipynb new file mode 100644 index 0000000..80e941b --- /dev/null +++ b/Notebooks/Beanies_processing.ipynb @@ -0,0 +1,1726 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Imports" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from skimage import measure\n", + "from scipy.interpolate import RegularGridInterpolator\n", + "from scipy.signal import savgol_filter\n", + "from scipy.optimize import curve_fit\n", + "import open3d as o3d\n", + "import matplotlib.animation as animation\n", + "from IPython.display import HTML\n", + "import matplotlib\n", + "from ipywidgets import interactive, fixed, IntSlider\n", + "import math\n", + "from skimage.morphology import skeletonize\n", + "from matplotlib.widgets import RectangleSelector\n", + "import nibabel as nib\n", + "from matplotlib.ticker import MaxNLocator\n", + "from plyfile import PlyData" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Helper Functions" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "# Function to marching cube meshes\n", + "def plotMesh(ax, mesh, ax_x, ax_y, ax_z, azim, elev):\n", + " \n", + " ax.add_collection3d(mesh)\n", + " ax.set_xlabel(\"x\")\n", + " ax.set_ylabel(\"y\")\n", + " ax.set_zlabel(\"z\")\n", + " ax.set_xlim(ax_x[0], ax_x[1]) \n", + " ax.set_ylim(ax_y[0], ax_y[1]) \n", + " ax.set_zlim(ax_z[0], ax_z[1]) \n", + " ax.set_aspect('equal')\n", + " ax.azim = azim\n", + " ax.elev = elev\n", + "\n", + "# Define the 3D polynomial function of 4th degree\n", + "def poly_3d(xy, a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13, a14, a15):\n", + " x, y = xy\n", + " return (a0 + a1*x + a2*y + a3*x**2 + a4*y**2 + a5*x*y + a6*x**3 + a7*y**3 + a8*x**2*y + a9*x*y**2 +\n", + " a10*x**4 + a11*y**4 + a12*x**3*y + a13*x**2*y**2 + a14*x*y**3)\n", + "\n", + "# Fits polynomial to set of datapoints and displays normalvectors (plot)\n", + "def fit_polynomial_surface(points, subset=10, normal_direction=0, num_points_quiver=10, num_points_surface=50, scale_factor=35):\n", + " '''\n", + " num_points_quiver = 10 # Adjust the number of points as needed for the quivers\n", + " num_points_surface = 50 # Adjust the number of points as needed for the surface\n", + " scale_factor = 35 # Adjust the scale of the normal vectors\n", + " normal_direction = 1 # Change this to 1 (or 0) to reverse the direction of the normals\n", + " WANT NORMALS FACING INTO BEANIE SKULL\n", + " \n", + " extra_distance: The extra distance added around the points for sampling\n", + " step_size: The step size between each sampling along the normal vectors\n", + " steps: The amount of steps done for sampling along normal vectors\n", + " grid_size_x: Size of the grid of sample points along the x-axis (can be left out and calculated automatically)\n", + " grid_size_y: Size of the grid of sample points along the y-axis (can be left out and calculated automatically)\n", + " '''\n", + "\n", + " # Extract x, y, z coordinates from the points (These are subsets defined earlier!)\n", + " x_num = points[::subset, 0]\n", + " y_num = points[::subset, 1]\n", + " z_num = points[::subset, 2]\n", + "\n", + " # Initial guess for the parameters\n", + " p0 = np.ones(16) # Needs to match the amount of a's in poly_3d (helper function)\n", + "\n", + " # Use curve_fit to fit the function to the data\n", + " popt, pcov = curve_fit(poly_3d, (x_num, y_num), z_num, p0)\n", + "\n", + " # Create a grid of x, y values for the surface\n", + " x_grid_surface, y_grid_surface = np.meshgrid(np.linspace(min(x_num), max(x_num), num_points_surface), np.linspace(min(y_num), max(y_num), num_points_surface))\n", + "\n", + " # Compute the corresponding z values for the surface\n", + " z_grid_surface = poly_3d((x_grid_surface, y_grid_surface), *popt)\n", + "\n", + " # Create a grid of x, y values for the quivers\n", + " x_grid_quiver, y_grid_quiver = np.meshgrid(np.linspace(min(x_num), max(x_num), num_points_quiver), np.linspace(min(y_num), max(y_num), num_points_quiver))\n", + "\n", + " # Compute the corresponding z values for the quivers\n", + " z_grid_quiver = poly_3d((x_grid_quiver, y_grid_quiver), *popt)\n", + "\n", + " # Define the symbols\n", + " a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13, a14, a15 = sp.symbols('a0 a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 a11 a12 a13 a14 a15')\n", + " x, y = sp.symbols('x y')\n", + "\n", + " # Define the polynomial function\n", + " poly_expr = (a0 + a1*x + a2*y + a3*x**2 + a4*y**2 + a5*x*y + a6*x**3 + a7*y**3 + a8*x**2*y + a9*x*y**2 +\n", + " a10*x**4 + a11*y**4 + a12*x**3*y + a13*x**2*y**2 + a14*x*y**3)\n", + "\n", + " # Differentiate with respect to x and y\n", + " dz_dx = sp.diff(poly_expr, x)\n", + " dz_dy = sp.diff(poly_expr, y)\n", + "\n", + " # Create lambdified functions for evaluating the derivatives\n", + " dz_dx_func = sp.lambdify((x, y, a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13, a14, a15), dz_dx)\n", + " dz_dy_func = sp.lambdify((x, y, a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13, a14, a15), dz_dy)\n", + "\n", + " # Evaluate the derivatives at each point on the quiver surface\n", + " dz_dx_vals = dz_dx_func(x_grid_quiver, y_grid_quiver, *popt)\n", + " dz_dy_vals = dz_dy_func(x_grid_quiver, y_grid_quiver, *popt)\n", + "\n", + " # Compute the normal vectors\n", + " normals = np.dstack((-dz_dx_vals, -dz_dy_vals, np.ones_like(dz_dx_vals)))\n", + "\n", + " # To change direction of normals\n", + " normals *= (-1)**normal_direction \n", + "\n", + " # Normalize the normals\n", + " normals /= np.linalg.norm(normals, axis=2, keepdims=True)\n", + "\n", + " # Create a 3D plot\n", + " fig = plt.figure()\n", + " ax = fig.add_subplot(111, projection='3d')\n", + "\n", + " # Plot the original points\n", + " ax.scatter(x_num, y_num, z_num, color='b')\n", + "\n", + " # Plot the fitted surface using the surface grid\n", + " ax.plot_surface(x_grid_surface, y_grid_surface, z_grid_surface, color='r', alpha=0.5)\n", + "\n", + " # Plot the normal vectors using the quiver grid\n", + " for i in range(len(x_grid_quiver)):\n", + " for j in range(len(y_grid_quiver)):\n", + " ax.quiver(x_grid_quiver[i, j], y_grid_quiver[i, j], z_grid_quiver[i, j], normals[i, j, 0], normals[i, j, 1], normals[i, j, 2],\n", + " length=scale_factor, linewidth=5, color='g')\n", + "\n", + " ax.set_aspect('equal')\n", + " plt.show()\n", + " \n", + " return popt, dz_dx_func, dz_dy_func\n", + "\n", + "# Generates PCA-coordinates to sample \n", + "def generate_sample_points(points, popt, dz_dx_func, dz_dy_func, normal_direction=0, subset=10, extra_distance=15, step_size=0.1, min_steps=10, grid_size_x=None, grid_size_y=None):\n", + " '''\n", + " subset: amount of subsamples used\n", + " extra_distance: The extra distance added around the points for sampling\n", + " step_size: The step size between each sampling along the normal vectors\n", + " steps: The amount of steps done for sampling along normal vectors\n", + " grid_size_x: Size of the grid of sample points along the x-axis (can be left out and calculated automatically)\n", + " grid_size_y: Size of the grid of sample points along the y-axis (can be left out and calculated automatically)\n", + "\n", + " STANDARD FOR AGAMODON:\n", + " extra_distance = 15\n", + " step_size = 0.1\n", + " min_steps: The number of steps in the least sampled direction! (The other direction is sampled 5x as much)\n", + " grid_size_x = 500\n", + " grid_size_y = 200\n", + " '''\n", + " # Extract x, y, z coordinates from the points (These are subsets defined earlier!)\n", + " x_num = points[::subset, 0]\n", + " y_num = points[::subset, 1]\n", + " z_num = points[::subset, 2]\n", + "\n", + " # For calculating the span of the points to generate meshgrid\n", + " min_values = np.min(points, axis=0)\n", + " max_values = np.max(points, axis=0)\n", + "\n", + " # If user has not chosen, make it the size of the spans\n", + " if grid_size_x is None:\n", + " grid_size_x = int(max_values[0] - min_values[0])\n", + " if grid_size_y is None:\n", + " grid_size_y = int(max_values[1] - min_values[1])\n", + " \n", + " # Step 1: Create a grid in the x-y-plane\n", + " x_grid, y_grid = np.meshgrid(np.linspace(min(x_num) - extra_distance, max(x_num) + extra_distance, grid_size_x), np.linspace(min(y_num) - extra_distance, max(y_num) + extra_distance, grid_size_y))\n", + "\n", + " # Step 2: Compute the height of the fitted polynomial surface at each point\n", + " z_grid = poly_3d((x_grid, y_grid), *popt)\n", + "\n", + " # Step 3: Compute the normal vector at each point on the grid\n", + " dz_dx_vals = dz_dx_func(x_grid, y_grid, *popt) # defined priorly. \n", + " dz_dy_vals = dz_dy_func(x_grid, y_grid, *popt)\n", + " normals = np.dstack((-dz_dx_vals, -dz_dy_vals, np.ones_like(dz_dx_vals)))\n", + "\n", + " normals *= (-1)**normal_direction # To change direction of normals\n", + "\n", + " # Step 4: Normalize the normals for better visualization\n", + " normals /= np.linalg.norm(normals, axis=2, keepdims=True)\n", + "\n", + " # Step 5: Go `steps` steps in the direction of the normal vector and `steps` steps in the opposite direction to create the layers\n", + " # Modify the number of steps in each direction\n", + " steps_positive = 5 * min_steps\n", + " steps_negative = min_steps\n", + " grid_3d = np.zeros((grid_size_y, grid_size_x, steps_positive + steps_negative + 1, 3)) # 3D array to store the x, y, and z coordinates of the points (each index in the 3D matrix contains a 3D-coordinate in PCA-space)\n", + "\n", + " # Go steps_negative steps in the opposite direction\n", + " for i in range(-steps_negative, 0):\n", + " displacement = i * step_size * normals # displacement in x, y, and z directions\n", + " grid_3d[:, :, i+steps_negative, 0] = x_grid + displacement[:, :, 0] # x coordinate\n", + " grid_3d[:, :, i+steps_negative, 1] = y_grid + displacement[:, :, 1] # y coordinate\n", + " grid_3d[:, :, i+steps_negative, 2] = z_grid + displacement[:, :, 2] # z coordinate\n", + "\n", + " # Go steps_positive steps in the direction of the normal vector\n", + " for i in range(steps_positive):\n", + " displacement = i * step_size * normals # displacement in x, y, and z directions\n", + " grid_3d[:, :, i+steps_negative+1, 0] = x_grid + displacement[:, :, 0] # x coordinate\n", + " grid_3d[:, :, i+steps_negative+1, 1] = y_grid + displacement[:, :, 1] # y coordinate\n", + " grid_3d[:, :, i+steps_negative+1, 2] = z_grid + displacement[:, :, 2] # z coordinate\n", + "\n", + " # Flatten the 3D grid\n", + " points_flat = grid_3d.reshape(-1, 3)\n", + "\n", + " # Plot a small subset of the found points\n", + " fig = plt.figure()\n", + " ax = fig.add_subplot(111, projection='3d')\n", + " ax.scatter(points_flat[::2000, 0], points_flat[::2000, 1], points_flat[::2000, 2])\n", + " plt.show()\n", + "\n", + " return points_flat, grid_3d.shape[0:3]\n", + "\n", + "# Interpolate sample points in original volume file\n", + "def sample_in_original_volume(points_original, original_nii_path, grid_shape, intMethod='linear', expVal=0.0):\n", + " \n", + " # Load in original volume\n", + " original = nib.load(original_nii_path)\n", + " imgDim = original.header['dim'][1:4]\n", + "\n", + " # Set-up interpolators for moving image\n", + " x = np.arange(start=0, stop=imgDim[0], step=1)\n", + " y = np.arange(start=0, stop=imgDim[1], step=1)\n", + " z = np.arange(start=0, stop=imgDim[2], step=1)\n", + " F_moving = RegularGridInterpolator((x, y, z), original.get_fdata().astype('float16'), method=intMethod, bounds_error=False, fill_value=expVal)\n", + "\n", + " # Evaluate transformed grid points in the moving image\n", + " fVal = F_moving(points_original)\n", + " # Reshape to voxel grid\n", + " volQ = np.reshape(fVal,newshape=grid_shape).astype('float32')\n", + "\n", + " return volQ\n", + "\n", + "# Animate moving through the layers of volume\n", + "def ani_through_volume(volume, indices=None):\n", + " '''\n", + " indices: List [a,b] of layer numbers you want animated\n", + " '''\n", + " def update_img(z):\n", + " ax.clear()\n", + " ax.imshow(volume[:,:,z])\n", + " ax.set_title('Layer ' + str(z))\n", + "\n", + " # User can specify indices they want to see\n", + " if indices != None:\n", + " a = indices[0]\n", + " b = indices[1]\n", + " else:\n", + " a = 0\n", + " b = volume.shape[2] - 1\n", + "\n", + " matplotlib.rcParams['animation.embed_limit'] = 2**128\n", + "\n", + " fig, ax = plt.subplots()\n", + " ani = animation.FuncAnimation(fig, update_img, frames=range(a,b+1), interval=200)\n", + " \n", + " return HTML(ani.to_jshtml())\n", + "\n", + "# Create interactive plot for slicing\n", + "def interactive_plot_slicing(volume):\n", + " def get_slice(vol, x_start, x_end, y_start, y_end, z):\n", + " slice = vol[x_start:x_end, y_start:y_end, z]\n", + " plt.imshow(slice)\n", + " plt.show()\n", + "\n", + " x_start_slider = IntSlider(min=0, max=volume.shape[0]-1, value=0)\n", + " x_end_slider = IntSlider(min=0, max=volume.shape[0]-1, value=volume.shape[0]-1)\n", + " y_start_slider = IntSlider(min=0, max=volume.shape[1]-1, value=0)\n", + " y_end_slider = IntSlider(min=0, max=volume.shape[1]-1, value=volume.shape[1]-1)\n", + " z_slider = IntSlider(min=0, max=volume.shape[2]-1, value=0)\n", + "\n", + " def update_x_end_range(*args):\n", + " x_end_slider.min = x_start_slider.value\n", + " x_start_slider.observe(update_x_end_range, 'value')\n", + "\n", + " def update_y_end_range(*args):\n", + " y_end_slider.min = y_start_slider.value\n", + " y_start_slider.observe(update_y_end_range, 'value')\n", + "\n", + " interactive_plot = interactive(get_slice, \n", + " vol=fixed(volume), \n", + " x_start=x_start_slider, \n", + " x_end=x_end_slider, \n", + " y_start=y_start_slider, \n", + " y_end=y_end_slider, \n", + " z=z_slider)\n", + " return interactive_plot\n", + "\n", + "# Allows user to remove unwanted sutures\n", + "def interactive_process_volume(volume, interactive_plot, z_vals):\n", + " # Define the slicing parameters\n", + " x_start = interactive_plot.kwargs['x_start']\n", + " x_end = interactive_plot.kwargs['x_end']\n", + " y_start = interactive_plot.kwargs['y_start']\n", + " y_end = interactive_plot.kwargs['y_end']\n", + " \n", + " min_z, max_z = z_vals[0], z_vals[1]\n", + "\n", + " # Crop the volume based on the provided slices\n", + " cropped_volume = volume[x_start:x_end, y_start:y_end, min_z:max_z+1]\n", + "\n", + " # Calculate the median layer index\n", + " median_z = cropped_volume.shape[2] // 2\n", + " median_layer = cropped_volume[:, :, median_z]\n", + "\n", + " # Create a new volume to apply changes\n", + " new_volume = cropped_volume.copy()\n", + "\n", + " fig, ax = plt.subplots()\n", + " ax.imshow(median_layer, cmap='gray')\n", + "\n", + " def line_select_callback(eclick, erelease):\n", + " nonlocal new_volume\n", + " x1, y1 = int(eclick.xdata), int(eclick.ydata)\n", + " x2, y2 = int(erelease.xdata), int(erelease.ydata)\n", + " if x1 > x2:\n", + " x1, x2 = x2, x1\n", + " if y1 > y2:\n", + " y1, y2 = y2, y1\n", + "\n", + " # Find the maximum value within the selected area\n", + " max_value = np.max(new_volume[y1:y2+1, x1:x2+1, :])\n", + "\n", + " # Assign the maximum value to all pixels within the selected area for each slice\n", + " for z in range(new_volume.shape[2]):\n", + " new_volume[y1:y2+1, x1:x2+1, z] = max_value\n", + "\n", + " # Update the displayed image in real-time\n", + " ax.clear()\n", + " ax.imshow(new_volume[:, :, median_z], cmap='gray')\n", + " plt.draw()\n", + "\n", + " rs = RectangleSelector(ax, line_select_callback,\n", + " useblit=True,\n", + " button=[1], # Only use left button\n", + " minspanx=5, minspany=5,\n", + " spancoords='pixels',\n", + " interactive=True)\n", + "\n", + " plt.show(block=True) # Block until plot window is closed\n", + "\n", + " # Return the new_volume after the plot is closed\n", + " return new_volume\n", + "\n", + "# Interactive plot for deciding threshold value\n", + "def interactive_plot_threshold(volume):\n", + "\n", + " # Define the function to generate the plot\n", + " def plot_threshold(threshold):\n", + " slice = volume[:,:, volume.shape[2]//2]\n", + " slice_bin = np.where(slice >= threshold, 1, 0)\n", + " slice_bin= 1 - slice_bin # Make suture the \"foreground\"\n", + " plt.imshow(slice_bin)\n", + " plt.show()\n", + "\n", + " # Calculate the minimum value for the slider\n", + " min_val = (volume.max()+volume.min())//2\n", + " # Round up to the nearest hundred\n", + " min_val = math.ceil(min_val / 100.0) * 100\n", + "\n", + " # Create a slider for the threshold value\n", + " threshold_slider = IntSlider(min=min_val, max=volume.max(), step=100, value=min_val)\n", + "\n", + " # Create the interactive plot\n", + " interactive_plot_new = interactive(plot_threshold, threshold=threshold_slider)\n", + "\n", + " return interactive_plot_new\n", + "\n", + "# Extracts the suture depending on your former choices\n", + "def extract_suture(volume, points_orig, interactive_plot_slice, interactive_plot_thresh, z_vals, grid_shape):\n", + " '''\n", + " Input volume is sliced\n", + " '''\n", + " min_z = z_vals[0]\n", + " max_z = z_vals[1]\n", + " num_indices = max_z - min_z\n", + "\n", + " # Extrac chosen indices\n", + " x_start = interactive_plot_slice.kwargs['x_start']\n", + " x_end = interactive_plot_slice.kwargs['x_end']\n", + " y_start = interactive_plot_slice.kwargs['y_start']\n", + " y_end = interactive_plot_slice.kwargs['y_end']\n", + "\n", + " # Extract chosen threshold\n", + " threshold_value = interactive_plot_thresh.kwargs['threshold']\n", + " \n", + " # Generate a slice\n", + " slice = volume[:,:, volume.shape[2]//2]\n", + "\n", + " # Generate empty array to store suture points in\n", + " suture = np.zeros((slice.shape[0], slice.shape[1], volume.shape[2]))\n", + " # For animation\n", + " binary_slices = np.zeros((slice.shape[0], slice.shape[1], num_indices))\n", + "\n", + " # Two first indicies were bad in segmentation\n", + " for index in range(num_indices):\n", + " slice = volume[:,:, index]\n", + " slice_bin = np.where(slice >= threshold_value, 1, 0)\n", + " slice_bin= 1 - slice_bin # Make suture the \"foreground\"\n", + "\n", + " # Label the connected components in the binary image\n", + " labels = measure.label(slice_bin)\n", + "\n", + " # Now, you can analyze each individual blob\n", + " properties = measure.regionprops(labels)\n", + "\n", + " # Sort the regions by area_bbox in descending order\n", + " properties_sorted = sorted(properties, key=lambda x: x.bbox_area, reverse=True)\n", + "\n", + " # Get the region with the largest area_bbox\n", + " largest_bbox_region = properties_sorted[0]\n", + "\n", + " # Get the coordinates of the pixels in the blob with the largest area_bbox\n", + " coords_largest_bbox = largest_bbox_region.coords\n", + " \n", + " # Extract x and y coordinates\n", + " x_coords = coords_largest_bbox[:, 1]\n", + " y_coords = coords_largest_bbox[:, 0]\n", + "\n", + " suture[y_coords, x_coords, index] = 1\n", + " binary_slices[:,:, index] = slice_bin\n", + "\n", + " ### ANIMATION\n", + " # Create a figure with two subplots\n", + " fig, axs = plt.subplots(1, 2, figsize=(8,3))\n", + "\n", + " # Initialize the images\n", + " im1 = axs[0].imshow(suture[:, :, 0], animated=True)\n", + " im2 = axs[1].imshow(volume[:,:, 0], animated=True)\n", + "\n", + " # Initialize the text\n", + " txt = axs[0].text(0.5, 1.01, f'Layer: {0}', transform=axs[0].transAxes)\n", + "\n", + " # Update function for the animation\n", + " def updatefig(i):\n", + " im1.set_array(suture[:, :, i])\n", + " im2.set_array(volume[:,:, i])\n", + " txt.set_text(f'Layer: {i}')\n", + " return im1, im2, txt\n", + "\n", + " # Create the animation\n", + " ani = animation.FuncAnimation(fig, updatefig, frames=range(num_indices), interval=200, blit=True)\n", + " \n", + " ### GET ORIGINAL COORDINATES\n", + " new_shape = grid_shape + (3,)\n", + " coords_orig = points_orig.reshape(new_shape)\n", + "\n", + " # cooresponding original coordinates to sliced area\n", + " coords_orig_sliced = coords_orig[x_start:x_end, y_start:y_end, min_z:min_z+num_indices]\n", + " \n", + " ### PLOT ORIGINAL COORDINATES OF SLICED AREA\n", + " # Get the indices where suture is 1\n", + " indices = np.where(suture == 1)\n", + "\n", + " # Get the corresponding points in coords_orig_sliced\n", + " points = coords_orig_sliced[indices]\n", + "\n", + " # Create a 3D scatter plot\n", + " fig = plt.figure()\n", + " ax = fig.add_subplot(111, projection='3d')\n", + " ax.scatter(points[::10, 0], points[::10, 1], points[::10, 2])\n", + " ax.set_xlabel('x')\n", + " ax.set_ylabel('y')\n", + " ax.set_zlabel('z')\n", + " plt.show()\n", + "\n", + " return ani, points, suture, coords_orig_sliced\n", + "\n", + "# Curve analysis\n", + "#############################################################\n", + "\n", + "def find_radius_3d(point1, point2, point3):\n", + " \"\"\"\n", + " Returns the center and radius of the circumsphere of a triangle in 3D space.\n", + " The circumsphere of a triangle is the sphere that passes through all three vertices of the triangle.\n", + " The function uses the formula from this link: https://gamedev.stackexchange.com/questions/60630/how-do-i-find-the-circumcenter-of-a-triangle-in-3d\n", + " \"\"\"\n", + "\n", + " # Convert the input points to numpy arrays for easier manipulation\n", + " p1 = np.array(point1)\n", + " p2 = np.array(point2)\n", + " p3 = np.array(point3)\n", + "\n", + " # Calculate the vectors from point1 to point2 and point1 to point3\n", + " p12 = p2 - p1 # Vector from point1 to point2\n", + " p13 = p3 - p1 # Vector from point1 to point3\n", + "\n", + " # Calculate the cross product of p12 and p13\n", + " p12_X_p13 = np.cross(p12, p13) # Cross product of vectors p12 and p13\n", + "\n", + " # Calculate the vector from point1 to the circumsphere center\n", + " toCircumsphereCenter = ((np.cross(p12_X_p13, p12) * np.dot(p13, p13)) + (np.cross(p13, p12_X_p13) * np.dot(p12, p12))) / (2 * np.dot(p12_X_p13, p12_X_p13))\n", + "\n", + " # Calculate the radius of the circumsphere\n", + " circumsphereRadius = np.linalg.norm(toCircumsphereCenter) # The radius is the length of the vector to the circumsphere center\n", + "\n", + " # Calculate the coordinates of the circumsphere center\n", + " ccs = point1 + toCircumsphereCenter # The center is point1 plus the vector to the circumsphere center\n", + "\n", + " # Return the center and radius of the circumsphere\n", + " return ccs, circumsphereRadius\n", + "\n", + "def my_scatter(points, camera_view=(30, -60)):\n", + " # Create a new figure\n", + " fig = plt.figure()\n", + "\n", + " # Create a 3D plot\n", + " ax = fig.add_subplot(111, projection='3d')\n", + "\n", + " # Plot the mean points\n", + " ax.scatter(points[:, 0], points[:, 1], points[:, 2]) # points\n", + "\n", + " # Set the labels with increased size\n", + " ax.set_xlabel('X', labelpad=20, fontsize=14)\n", + " ax.set_ylabel('Y', labelpad=20, fontsize=14)\n", + " ax.set_zlabel('Z', labelpad=20, fontsize=14)\n", + "\n", + " # Increase the size of the tick labels\n", + " ax.tick_params(axis='both', which='major', labelsize=14)\n", + "\n", + " # Set the camera view\n", + " ax.view_init(*camera_view)\n", + "\n", + " #print(\"Number of points: \", points.shape[0])\n", + " plt.show()\n", + " return None\n", + "\n", + "def load_data_ply(file):\n", + " plydata = PlyData.read(file)\n", + "\n", + " # Extract the vertex data\n", + " vertex_data = plydata['vertex'].data\n", + "\n", + " # Extract x, y, z coordinates\n", + " x = vertex_data['x']\n", + " y = vertex_data['y']\n", + " z = vertex_data['z']\n", + "\n", + " vertex_data = np.column_stack([x,y,z])\n", + " return vertex_data\n", + "\n", + "def load_data_npy(file):\n", + " # Load the data from the .npy file\n", + " vertex_data = np.load(file)\n", + " return vertex_data\n", + "\n", + "def mean_points_func(vertex_data, plot = False, N=4):\n", + " # This is from ply data, we chose resolution 0, so its a square (4 points)\n", + "\n", + " # Initialize an empty list to store the mean points\n", + " mean_points = []\n", + "\n", + " # Loop over the vertex_data array with a step of 4\n", + " for i in range(0, vertex_data.shape[0] - 3, N):\n", + " # Select the current 4 points\n", + " current_points = vertex_data[i:i+N]\n", + " \n", + " # Compute the mean of the current 4 points\n", + " mean_point = current_points.mean(axis=0)\n", + " \n", + " # Append the mean point to the list\n", + " mean_points.append(mean_point)\n", + "\n", + " # Convert the list to a numpy array\n", + " mean_points = np.array(mean_points)\n", + " if plot == True:\n", + " my_scatter(mean_points)\n", + " return mean_points\n", + "\n", + "def my_sort_points(points, threshold):\n", + " \"\"\"\n", + " This function sorts points based on their Euclidean distance. Starting from the point with the lowest x value,\n", + " it finds the next closest point and adds it to the sorted list. This process continues until no points are left,\n", + " or the distance to the next closest point is greater than a specified threshold.\n", + "\n", + " Parameters:\n", + " points (numpy.ndarray): The points to be sorted. Each point is an array of coordinates.\n", + " threshold (Int): The maximum allowed distance to the next point.\n", + "\n", + " Returns:\n", + " numpy.ndarray: The sorted points.\n", + " \"\"\"\n", + "\n", + " # Find the point with the lowest x value\n", + " start_point = min(points, key=lambda point: point[0])\n", + "\n", + " # Initialize the sorted points list with the start point\n", + " sorted_points = [tuple(start_point)]\n", + "\n", + " # Create a set of the remaining points\n", + " remaining_points = set(map(tuple, points))\n", + " # Remove the start point from the remaining points\n", + " remaining_points.remove(tuple(start_point))\n", + "\n", + " # Continue until no points are left\n", + " while remaining_points:\n", + " # Get the last point added to the sorted list\n", + " current_point = sorted_points[-1]\n", + "\n", + " # Calculate the Euclidean distance from the current point to each remaining point\n", + " distances = {point: np.linalg.norm(np.array(point) - np.array(current_point)) for point in remaining_points}\n", + "\n", + " # Find the point with the minimum distance to the current point\n", + " next_point, min_distance = min(distances.items(), key=lambda item: item[1])\n", + "\n", + " # If the minimum distance is greater than the threshold, stop sorting\n", + " if min_distance > threshold:\n", + " break\n", + "\n", + " # Remove the next point from the remaining points\n", + " remaining_points.remove(next_point)\n", + " # Add the next point to the sorted list\n", + " sorted_points.append(next_point)\n", + "\n", + " # Return the sorted points as a numpy array\n", + " return np.array(sorted_points)\n", + "\n", + "def my_sort_points_max_z(points, threshold):\n", + " \"\"\"\n", + " This function sorts points based on their Euclidean distance. Starting from the point with the lowest x value,\n", + " it finds the next closest point and adds it to the sorted list. This process continues until no points are left,\n", + " or the distance to the next closest point is greater than a specified threshold.\n", + "\n", + " Parameters:\n", + " points (numpy.ndarray): The points to be sorted. Each point is an array of coordinates.\n", + " threshold (Int): The maximum allowed distance to the next point.\n", + "\n", + " Returns:\n", + " numpy.ndarray: The sorted points.\n", + " \"\"\"\n", + "\n", + " # Find the point with the lowest x value\n", + " start_point = max(points, key=lambda point: point[2])\n", + "\n", + " # Initialize the sorted points list with the start point\n", + " sorted_points = [tuple(start_point)]\n", + "\n", + " # Create a set of the remaining points\n", + " remaining_points = set(map(tuple, points))\n", + " # Remove the start point from the remaining points\n", + " remaining_points.remove(tuple(start_point))\n", + "\n", + " # Continue until no points are left\n", + " while remaining_points:\n", + " # Get the last point added to the sorted list\n", + " current_point = sorted_points[-1]\n", + "\n", + " # Calculate the Euclidean distance from the current point to each remaining point\n", + " distances = {point: np.linalg.norm(np.array(point) - np.array(current_point)) for point in remaining_points}\n", + "\n", + " # Find the point with the minimum distance to the current point\n", + " next_point, min_distance = min(distances.items(), key=lambda item: item[1])\n", + "\n", + " # If the minimum distance is greater than the threshold, stop sorting\n", + " if min_distance > threshold:\n", + " break\n", + "\n", + " # Remove the next point from the remaining points\n", + " remaining_points.remove(next_point)\n", + " # Add the next point to the sorted list\n", + " sorted_points.append(next_point)\n", + "\n", + " # Return the sorted points as a numpy array\n", + " return np.array(sorted_points)\n", + "\n", + "def my_sort_points_max_x(points, threshold):\n", + " \"\"\"\n", + " This function sorts points based on their Euclidean distance. Starting from the point with the lowest x value,\n", + " it finds the next closest point and adds it to the sorted list. This process continues until no points are left,\n", + " or the distance to the next closest point is greater than a specified threshold.\n", + "\n", + " Parameters:\n", + " points (numpy.ndarray): The points to be sorted. Each point is an array of coordinates.\n", + " threshold (Int): The maximum allowed distance to the next point.\n", + "\n", + " Returns:\n", + " numpy.ndarray: The sorted points.\n", + " \"\"\"\n", + "\n", + " # Find the point with the lowest x value\n", + " start_point = max(points, key=lambda point: point[0])\n", + "\n", + " # Initialize the sorted points list with the start point\n", + " sorted_points = [tuple(start_point)]\n", + "\n", + " # Create a set of the remaining points\n", + " remaining_points = set(map(tuple, points))\n", + " # Remove the start point from the remaining points\n", + " remaining_points.remove(tuple(start_point))\n", + "\n", + " # Continue until no points are left\n", + " while remaining_points:\n", + " # Get the last point added to the sorted list\n", + " current_point = sorted_points[-1]\n", + "\n", + " # Calculate the Euclidean distance from the current point to each remaining point\n", + " distances = {point: np.linalg.norm(np.array(point) - np.array(current_point)) for point in remaining_points}\n", + "\n", + " # Find the point with the minimum distance to the current point\n", + " next_point, min_distance = min(distances.items(), key=lambda item: item[1])\n", + "\n", + " # If the minimum distance is greater than the threshold, stop sorting\n", + " if min_distance > threshold:\n", + " break\n", + "\n", + " # Remove the next point from the remaining points\n", + " remaining_points.remove(next_point)\n", + " # Add the next point to the sorted list\n", + " sorted_points.append(next_point)\n", + "\n", + " # Return the sorted points as a numpy array\n", + " return np.array(sorted_points)\n", + "\n", + "def interpolate_3d_points(points, num_interpolated_points = 3, skip_points=5, camera_view=(30, -60)):\n", + " \"\"\"\n", + " This function interpolates additional points between each pair of consecutive points in a 3D space.\n", + " It then skips a specified number of points to reduce the density of points.\n", + "\n", + " Parameters:\n", + " points (numpy.ndarray): The original points in 3D space.\n", + " num_interpolated_points (int): The number of points to interpolate between each pair of original points.\n", + " skip_points (int): The number of points to skip in the final list of points.\n", + "\n", + " Returns:\n", + " numpy.ndarray: The interpolated points.\n", + " \"\"\"\n", + "\n", + " # Create a list to store the new points\n", + " interpolated_points = []\n", + " \n", + " # Iterate over each pair of consecutive points\n", + " for i in range(len(points) - 1):\n", + " # Get the start and end points of the current segment\n", + " start_point = points[i]\n", + " end_point = points[i + 1]\n", + " \n", + " # Append the start point of the current segment to the list\n", + " interpolated_points.append(start_point)\n", + " \n", + " # Compute the interpolated points for the current segment\n", + " for j in range(1, num_interpolated_points + 1):\n", + " # Compute the parameter for the linear interpolation\n", + " t = j / (num_interpolated_points + 1)\n", + " # Compute the interpolated point\n", + " interpolated_point = start_point * (1 - t) + end_point * t\n", + " # Append the interpolated point to the list\n", + " interpolated_points.append(interpolated_point)\n", + " \n", + " # Append the last point of the original array to the list\n", + " interpolated_points.append(points[-1])\n", + "\n", + " # Skip points to reduce the density\n", + " interpolated_points = interpolated_points[::skip_points]\n", + "\n", + " # Print the number of interpolated points\n", + " print(\"Number of interpolated points: \", len(interpolated_points))\n", + "\n", + " # Visualize the interpolated points\n", + " my_scatter(np.array(interpolated_points), camera_view=camera_view)\n", + " # Return the interpolated points as a numpy array\n", + " return np.array(interpolated_points)\n", + "\n", + "def smoothing_signal(points, window_length, polyorder, camera_view=(30, -60)):\n", + " \"\"\"\n", + " This function applies a Savitzky-Golay filter to smooth a signal represented by a set of points.\n", + " It then visualizes the smoothed points using a scatter plot.\n", + "\n", + " Parameters:\n", + " points (numpy.ndarray): The original points representing the signal.\n", + " window_length (int): The length of the filter window (i.e., the number of points used for the polynomial regression).\n", + " polyorder (int): The order of the polynomial used in the Savitzky-Golay filter.\n", + "\n", + " Returns:\n", + " numpy.ndarray: The smoothed points.\n", + " \"\"\"\n", + "\n", + " # Apply the Savitzky-Golay filter to the points\n", + " # The filter fits a polynomial of a specified order to a window of points and uses the polynomial to estimate the center point of the window\n", + " # This process is repeated for each point in the signal, resulting in a smoothed version of the signal\n", + " smooth_points = savgol_filter(points, window_length, polyorder, axis=0)\n", + "\n", + " # Visualize the points\n", + " my_scatter(smooth_points, camera_view)\n", + "\n", + " # Print the number of smoothed points\n", + " print(\"Number of smoothed points: \", len(smooth_points))\n", + "\n", + " # Return the smoothed points\n", + " return smooth_points\n", + "\n", + "def calculate_radii(points):\n", + " \"\"\"\n", + " This function calculates the centers and radii of the circles passing through each triplet of consecutive points.\n", + " It then appends a large radius at the start and end of the radii array to represent the end points of the line.\n", + "\n", + " Parameters:\n", + " points (numpy.ndarray): The points through which the circles pass.\n", + "\n", + " Returns:\n", + " tuple: A tuple containing two numpy arrays. The first array contains the centers of the circles, and the second array contains the radii of the circles.\n", + " \"\"\"\n", + "\n", + " # Initialize lists to store the centers and radii\n", + " centers = []\n", + " radii = []\n", + "\n", + " # Loop over the points\n", + " for i in range(1, len(points) - 1):\n", + " # Find the circle passing through the points points[i-1], points[i], points[i+1]\n", + " center, radius = find_radius_3d(points[i-1], points[i], points[i+1])\n", + " \n", + " # Append the center and radius to the lists\n", + " centers.append(center)\n", + " radii.append(radius)\n", + "\n", + " # Convert the lists to numpy arrays\n", + " centers = np.array(centers)\n", + " radii = np.array(radii)\n", + "\n", + " # Append a large number at the end and start of the radii array to represent the end points of the line\n", + " radii = np.append(radii, 10**5)\n", + " radii = np.insert(radii, 0, 10**5)\n", + "\n", + " # Print the shapes of the radii and points arrays\n", + " print(f\"Radii shape: {radii.shape}\\nequidistant_points shape: {points.shape}\")\n", + " # Print the minimum and maximum of the radii\n", + " print(f\"min of radii: {np.min(radii)}\\nmax of radii: {np.max(radii)}\")\n", + "\n", + " # Return the centers and radii\n", + " return centers, radii\n", + "\n", + "def plot_thresholded_points(smooth_points, radii, T=18, camera_view=(30, -60)):\n", + " # Sort the radii\n", + " radii_sort = np.sort(radii)\n", + "\n", + " threshold = radii_sort[T] # use the sorted list since the lower values should be the curves\n", + " # Create a mask for the radii below the threshold\n", + " mask = radii < threshold # use original array to get the correct index for points\n", + "\n", + " # Select the points corresponding to the radii below the threshold\n", + " selected_points = smooth_points[mask]\n", + "\n", + " # Create a new figure\n", + " fig = plt.figure()\n", + "\n", + " # Create a 3D plot\n", + " ax = fig.add_subplot(111, projection='3d')\n", + "\n", + " # Plot the selected points\n", + " ax.scatter(smooth_points[:, 0], smooth_points[:, 1], smooth_points[:, 2], color = \"b\",alpha=0.5) \n", + " ax.scatter(selected_points[:, 0], selected_points[:, 1], selected_points[:, 2], color = \"r\",alpha=1, s = 50 , marker='o')\n", + " # Set the labels with increased size\n", + " ax.set_xlabel('X', labelpad=20, fontsize=14)\n", + " ax.set_ylabel('Y', labelpad=20, fontsize=14)\n", + " ax.set_zlabel('Z', labelpad=20, fontsize=14)\n", + "\n", + " # Increase the size of the tick labels\n", + " ax.tick_params(axis='both', which='major', labelsize=14)\n", + "\n", + " # Set the camera view\n", + " ax.view_init(*camera_view)\n", + "\n", + " print(f\"Total number of points below threshold: {len(selected_points)}\\n\")\n", + " # Show the plot\n", + " plt.show()\n", + " # må lige lege med at finde ud af hvad man vil have af output\n", + " return selected_points\n", + "\n", + "def plot_turning_points(smooth_points, radii, T = 18, skip_indices = 4, camera_view=(30, -60)):\n", + "\n", + " radii_sort = np.sort(radii)\n", + " # Define the threshold\n", + " threshold = radii_sort[T]\n", + "\n", + " # Create a mask for the radii below the threshold\n", + " # This will create a boolean array where each element is True \n", + " # if the corresponding radius is below the threshold, and False otherwise\n", + " mask = radii < threshold\n", + "\n", + " # Select the points corresponding to the radii below the threshold\n", + " # This will create a new array containing only the points that correspond to the radii below the threshold\n", + " selected_points = smooth_points[mask]\n", + "\n", + " # Initialize an empty list to keep track of the indices of the points that have been plotted\n", + " plotted_indices = []\n", + "\n", + " # Create a new figure\n", + " fig = plt.figure()\n", + "\n", + " # Add a 3D subplot to the figure\n", + " ax = fig.add_subplot(111, projection='3d')\n", + "\n", + " # Iterate over the selected points and their indices in the original points array\n", + " for point, index in zip(selected_points, np.where(mask)[0]):\n", + " # Check if the current index is within a range of 4 indices of any previously plotted point\n", + " # If it is not, plot the point and add its index to the list of plotted indices\n", + " if not any(abs(index - plotted_index) <= skip_indices for plotted_index in plotted_indices):\n", + " ax.scatter(point[0], point[1], point[2], color = \"r\", alpha=1, s=100, marker='o', edgecolors='black')\n", + " plotted_indices.append(index)\n", + "\n", + " # Plot the original points in blue\n", + " ax.scatter(smooth_points[:, 0], smooth_points[:, 1], smooth_points[:, 2], color = \"b\", alpha=0.5)\n", + "\n", + " # Set the labels with increased size\n", + " ax.set_xlabel('X', labelpad=20, fontsize=14)\n", + " ax.set_ylabel('Y', labelpad=20, fontsize=14)\n", + " ax.set_zlabel('Z', labelpad=20, fontsize=14)\n", + "\n", + " # Increase the size of the tick labels\n", + " ax.tick_params(axis='both', which='major', labelsize=14)\n", + "\n", + " # Set the camera view\n", + " ax.view_init(*camera_view)\n", + "\n", + " plt.show()\n", + " print(f\"Number of turning points: {len(plotted_indices)}\\n\")\n", + " return plotted_indices\n", + "\n", + "def curve_characteristics(points):\n", + " start_p, end_p = points[0],points[-1]\n", + " distance = np.linalg.norm(start_p - end_p)\n", + " print(f\" The distance between the start and end point is: {distance} Voxels\")\n", + " # Initialize a variable to store the total length\n", + " total_length = 0\n", + "\n", + " # Loop over the points\n", + " for i in range(len(points) - 1):\n", + " # Compute the Euclidean distance between the current point and the next point\n", + " distance2 = np.linalg.norm(points[i] - points[i+1])\n", + " \n", + " # Add the distance to the total length\n", + " total_length += distance2\n", + "\n", + " print(f\" The total length of the line is: {total_length} Voxels\")\n", + " print(f\" The tortuosity of the line is: {total_length/distance}\")\n", + "\n", + " # Assuming smoothed_points is a numpy array of shape (n, 3)\n", + " min_vals = np.min(points, axis=0)\n", + " max_vals = np.max(points, axis=0)\n", + "\n", + " # Calculate the lengths of the sides of the bounding box\n", + " lengths = max_vals - min_vals\n", + "\n", + " # Calculate the volume of the bounding box\n", + " volume = np.prod(lengths)\n", + "\n", + " print(\" The volume of the bounding box is: \", volume, \"Voxels\")\n", + " return distance, total_length, total_length/distance, volume\n", + "\n", + "############## Interactive plots of different functions ####################\n", + "\n", + "def interactive_interp_points(mean_points):\n", + " # Create sliders\n", + " num_interpolated_points_slider = widgets.IntSlider(min=0, max=20, step=1, value=5, description='Interpolated Points:')\n", + " skip_points_slider = widgets.IntSlider(min=1, max=40, step=1, value=20)\n", + " elevation_slider = widgets.IntSlider(min=0, max=90, step=1, value=30, description='Elevation:')\n", + " azimuth_slider = widgets.IntSlider(min=-180, max=180, step=1, value=-60, description='Azimuth:')\n", + "\n", + "\n", + " # Function to update plot\n", + " def update_plot(num_interpolated_points, skip_points, elevation, azimuth):\n", + " interp_points = interpolate_3d_points(mean_points, num_interpolated_points, skip_points, camera_view=(elevation, azimuth))\n", + " return interp_points\n", + "\n", + " # Interactive widget\n", + " interactive_plot = widgets.interactive(update_plot, num_interpolated_points=num_interpolated_points_slider, skip_points=skip_points_slider, elevation=elevation_slider, azimuth=azimuth_slider)\n", + " return interactive_plot\n", + "\n", + "def interactive_smooth_points(interp_points):\n", + " # Create sliders\n", + " window_length_slider = widgets.IntSlider(min=1, max=50, step=1, value=10, description='Window_length:')\n", + " polyorder_slider = widgets.IntSlider(min=1, max=8, step=1, value=1, description='Polyorder:')\n", + " elevation_slider = widgets.IntSlider(min=0, max=90, step=1, value=30, description='Elevation:')\n", + " azimuth_slider = widgets.IntSlider(min=-180, max=180, step=1, value=-60, description='Azimuth:')\n", + "\n", + " def update_plot(window_length, polyorder, elevation, azimuth):\n", + " if window_length <= polyorder:\n", + " print(\"Error: window_length must be greater than polyorder.\")\n", + " return\n", + " smoothed_points = smoothing_signal(interp_points, window_length, polyorder,camera_view=(elevation, azimuth))\n", + " \n", + " return smoothed_points\n", + " \n", + " # Interactive widget\n", + " interactive_smooth = widgets.interactive(update_plot, window_length=window_length_slider, polyorder=polyorder_slider, elevation=elevation_slider, azimuth=azimuth_slider)\n", + " return interactive_smooth\n", + "\n", + "def interactive_thresholed_points(smoothed_points, radii):\n", + " # Create sliders\n", + " T_slider = widgets.IntSlider(min=0, max=100, step=1, value=55, description='T:')\n", + " elevation_slider = widgets.IntSlider(min=0, max=90, step=1, value=30, description='Elevation:')\n", + " azimuth_slider = widgets.IntSlider(min=-180, max=180, step=1, value=-60, description='Azimuth:')\n", + "\n", + " # Function to update plot\n", + " def update_plot(T, elevation, azimuth):\n", + " selected_points = plot_thresholded_points(smoothed_points, radii, T,camera_view=(elevation, azimuth))\n", + " return selected_points\n", + "\n", + " # Interactive widget\n", + " interactive_selected_points = widgets.interactive(update_plot, T=T_slider, elevation=elevation_slider, azimuth=azimuth_slider)\n", + " return interactive_selected_points\n", + "\n", + "def interactive_turning_points(smoothed_points, radii):\n", + " # Create sliders\n", + " T_slider = widgets.IntSlider(min=0, max=100, step=1, value=60, description='T:')\n", + " skip_indices_slider = widgets.IntSlider(min=0, max=100, step=1, value=13, description='Skip Indices:')\n", + " elevation_slider = widgets.IntSlider(min=0, max=90, step=1, value=30, description='Elevation:')\n", + " azimuth_slider = widgets.IntSlider(min=-180, max=180, step=1, value=-60, description='Azimuth:')\n", + "\n", + " # Function to update plot\n", + " def update_plot(T, skip_indices, elevation, azimuth):\n", + " selected_turning_points = plot_turning_points(smoothed_points, radii, T, skip_indices, camera_view=(elevation, azimuth))\n", + " return selected_turning_points\n", + "\n", + " # Interactive widget\n", + " interactive_turning_points = widgets.interactive(update_plot, T=T_slider, skip_indices=skip_indices_slider, elevation=elevation_slider, azimuth=azimuth_slider)\n", + " return interactive_turning_points\n", + "\n", + "def interactive_sorting_points(vertex_data):\n", + "\n", + " # Create sliders\n", + " threshold_slider = widgets.FloatSlider(min=1, max=100, step=1, value=15, description='Threshold:')\n", + " elevation_slider = widgets.IntSlider(min=0, max=90, step=1, value=30, description='Elevation:')\n", + " azimuth_slider = widgets.IntSlider(min=-180, max=180, step=1, value=-60, description='Azimuth:')\n", + "\n", + " # Function to update plot\n", + " def update_plot(threshold, elevation, azimuth):\n", + " sorted_points = my_sort_points(vertex_data, threshold)\n", + " my_scatter(sorted_points,camera_view=(elevation, azimuth))\n", + " return sorted_points\n", + "\n", + " # Interactive widget\n", + " interactive_sorted_points = widgets.interactive(update_plot, threshold=threshold_slider, elevation=elevation_slider, azimuth=azimuth_slider)\n", + " return interactive_sorted_points\n", + "\n", + "def interactive_sorting_points_max_z(vertex_data):\n", + " # Create sliders\n", + " threshold_slider = widgets.FloatSlider(min=1, max=100, step=1, value=15, description='Threshold:')\n", + " elevation_slider = widgets.IntSlider(min=0, max=90, step=1, value=30, description='Elevation:')\n", + " azimuth_slider = widgets.IntSlider(min=-180, max=180, step=1, value=-60, description='Azimuth:')\n", + "\n", + " # Function to update plot\n", + " def update_plot(threshold, elevation, azimuth):\n", + " sorted_points = my_sort_points_max_z(vertex_data, threshold)\n", + " my_scatter(sorted_points,camera_view=(elevation, azimuth))\n", + " return sorted_points\n", + "\n", + " # Interactive widget\n", + " interactive_sorted_points = widgets.interactive(update_plot, threshold=threshold_slider, elevation=elevation_slider, azimuth=azimuth_slider)\n", + " return interactive_sorted_points\n", + "\n", + "def interactive_sorting_points_max_x(vertex_data):\n", + " # Create sliders\n", + " threshold_slider = widgets.FloatSlider(min=1, max=100, step=1, value=15, description='Threshold:')\n", + " elevation_slider = widgets.IntSlider(min=0, max=90, step=1, value=30, description='Elevation:')\n", + " azimuth_slider = widgets.IntSlider(min=-180, max=180, step=1, value=-60, description='Azimuth:')\n", + "\n", + " # Function to update plot\n", + " def update_plot(threshold, elevation, azimuth):\n", + " sorted_points = my_sort_points_max_x(vertex_data, threshold)\n", + " my_scatter(sorted_points,camera_view=(elevation, azimuth))\n", + " return sorted_points\n", + "\n", + " # Interactive widget\n", + " interactive_sorted_points = widgets.interactive(update_plot, threshold=threshold_slider, elevation=elevation_slider, azimuth=azimuth_slider)\n", + " return interactive_sorted_points\n", + "\n", + "def interactive_points(sorted_points):\n", + " fig = plt.figure()\n", + " ax = fig.add_subplot(111, projection='3d')\n", + "\n", + " scatter = ax.scatter([], [], [])\n", + "\n", + " def update(num):\n", + " ax.clear()\n", + " ax.scatter(sorted_points[:num+1, 0], sorted_points[:num+1, 1], sorted_points[:num+1, 2])\n", + " ax.set_xlim([sorted_points[:, 0].min(), sorted_points[:, 0].max()])\n", + " ax.set_ylim([sorted_points[:, 1].min(), sorted_points[:, 1].max()])\n", + " ax.set_zlim([sorted_points[:, 2].min(), sorted_points[:, 2].max()])\n", + "\n", + " ani = animation.FuncAnimation(fig, update, frames=len(sorted_points), interval=200)\n", + "\n", + " return HTML(ani.to_jshtml())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Unfolded volume generation" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[Open3D WARNING] Read PLY failed: unable to open file: YOUR_FILE_PATH.ply\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "c:\\Users\\aske_\\anaconda3\\envs\\Fagprojekt\\lib\\site-packages\\numpy\\core\\fromnumeric.py:3504: RuntimeWarning: Mean of empty slice.\n", + " return _methods._mean(a, axis=axis, dtype=dtype,\n", + "c:\\Users\\aske_\\anaconda3\\envs\\Fagprojekt\\lib\\site-packages\\numpy\\core\\_methods.py:121: RuntimeWarning: invalid value encountered in divide\n", + " ret = um.true_divide(\n", + "c:\\Users\\aske_\\anaconda3\\envs\\Fagprojekt\\lib\\site-packages\\numpy\\lib\\function_base.py:520: RuntimeWarning: Mean of empty slice.\n", + " avg = a.mean(axis, **keepdims_kw)\n", + "C:\\Users\\aske_\\AppData\\Local\\Temp\\ipykernel_25708\\3941942745.py:13: RuntimeWarning: Degrees of freedom <= 0 for slice\n", + " cov = np.cov(np.transpose(points)) # Covariance matrix\n", + "c:\\Users\\aske_\\anaconda3\\envs\\Fagprojekt\\lib\\site-packages\\numpy\\lib\\function_base.py:2748: RuntimeWarning: divide by zero encountered in divide\n", + " c *= np.true_divide(1, fact)\n", + "c:\\Users\\aske_\\anaconda3\\envs\\Fagprojekt\\lib\\site-packages\\numpy\\lib\\function_base.py:2748: RuntimeWarning: invalid value encountered in multiply\n", + " c *= np.true_divide(1, fact)\n" + ] + }, + { + "ename": "LinAlgError", + "evalue": "SVD did not converge", + "output_type": "error", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[1;31mLinAlgError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[1;32mIn[5], line 16\u001b[0m\n\u001b[0;32m 13\u001b[0m cov \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mcov(np\u001b[38;5;241m.\u001b[39mtranspose(points)) \u001b[38;5;66;03m# Covariance matrix \u001b[39;00m\n\u001b[0;32m 15\u001b[0m \u001b[38;5;66;03m# SVD or Eigendecomposition\u001b[39;00m\n\u001b[1;32m---> 16\u001b[0m U,S,V \u001b[38;5;241m=\u001b[39m \u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mlinalg\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43msvd\u001b[49m\u001b[43m(\u001b[49m\u001b[43mcov\u001b[49m\u001b[43m,\u001b[49m\u001b[43mfull_matrices\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mFalse\u001b[39;49;00m\u001b[43m)\u001b[49m\n\u001b[0;32m 18\u001b[0m \u001b[38;5;66;03m# Apply transform to rotated pointcloud\u001b[39;00m\n\u001b[0;32m 19\u001b[0m points_rot \u001b[38;5;241m=\u001b[39m (points \u001b[38;5;241m-\u001b[39m points_mu) \u001b[38;5;241m@\u001b[39m U\n", + "File \u001b[1;32mc:\\Users\\aske_\\anaconda3\\envs\\Fagprojekt\\lib\\site-packages\\numpy\\linalg\\linalg.py:1681\u001b[0m, in \u001b[0;36msvd\u001b[1;34m(a, full_matrices, compute_uv, hermitian)\u001b[0m\n\u001b[0;32m 1678\u001b[0m gufunc \u001b[38;5;241m=\u001b[39m _umath_linalg\u001b[38;5;241m.\u001b[39msvd_n_s\n\u001b[0;32m 1680\u001b[0m signature \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mD->DdD\u001b[39m\u001b[38;5;124m'\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m isComplexType(t) \u001b[38;5;28;01melse\u001b[39;00m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124md->ddd\u001b[39m\u001b[38;5;124m'\u001b[39m\n\u001b[1;32m-> 1681\u001b[0m u, s, vh \u001b[38;5;241m=\u001b[39m \u001b[43mgufunc\u001b[49m\u001b[43m(\u001b[49m\u001b[43ma\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43msignature\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43msignature\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mextobj\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mextobj\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 1682\u001b[0m u \u001b[38;5;241m=\u001b[39m u\u001b[38;5;241m.\u001b[39mastype(result_t, copy\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mFalse\u001b[39;00m)\n\u001b[0;32m 1683\u001b[0m s \u001b[38;5;241m=\u001b[39m s\u001b[38;5;241m.\u001b[39mastype(_realType(result_t), copy\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mFalse\u001b[39;00m)\n", + "File \u001b[1;32mc:\\Users\\aske_\\anaconda3\\envs\\Fagprojekt\\lib\\site-packages\\numpy\\linalg\\linalg.py:121\u001b[0m, in \u001b[0;36m_raise_linalgerror_svd_nonconvergence\u001b[1;34m(err, flag)\u001b[0m\n\u001b[0;32m 120\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_raise_linalgerror_svd_nonconvergence\u001b[39m(err, flag):\n\u001b[1;32m--> 121\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m LinAlgError(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mSVD did not converge\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n", + "\u001b[1;31mLinAlgError\u001b[0m: SVD did not converge" + ] + } + ], + "source": [ + "# Cropped mesh (as PLY-file)\n", + "input_file = 'YOUR_FILE.ply'\n", + "\n", + "pcd = o3d.io.read_point_cloud(input_file) \n", + "points = np.asarray(pcd.points)\n", + "\n", + "# Choose a subset of points and view point of display\n", + "subset = 10\n", + "viewPoint = [90, -45]\n", + "\n", + "# Principal Component Analysis (PCA) of the point cloud\n", + "points_mu = np.mean(points, axis = 0) # Center the point cloud\n", + "cov = np.cov(np.transpose(points)) # Covariance matrix \n", + "\n", + "# SVD or Eigendecomposition\n", + "U,S,V = np.linalg.svd(cov,full_matrices=False)\n", + "\n", + "# Apply transform to rotated pointcloud\n", + "points_rot = (points - points_mu) @ U\n", + "\n", + "## PLOT\n", + "# Plot the rotated pointcloud\n", + "fig = plt.figure(figsize=(8,4))\n", + "\n", + "# display the original point cloud with the principal axes\n", + "ax = fig.add_subplot(1,2,1, projection='3d')\n", + "ax.scatter(points[::subset,0], points[::subset,1], points[::subset,2], c='k', marker='.')\n", + "# Set the number of ticks on each axis\n", + "ax.xaxis.set_major_locator(MaxNLocator(nbins=5))\n", + "ax.yaxis.set_major_locator(MaxNLocator(nbins=5))\n", + "ax.zaxis.set_major_locator(MaxNLocator(nbins=5))\n", + "ax.set_xlabel(\"x\")\n", + "ax.set_ylabel(\"y\")\n", + "ax.set_zlabel(\"z\")\n", + "ax.set_aspect('equal')\n", + "ax.azim = 40\n", + "ax.elev = -25\n", + "ax.set_title('Original point cloud')\n", + "\n", + "# Add the principal axes to the plot using quiver\n", + "ax.quiver(points_mu[0], points_mu[1], points_mu[2], U[0,0], U[1,0], U[2,0], color='r', length=100, normalize=True, linewidth=2.5)\n", + "ax.quiver(points_mu[0], points_mu[1], points_mu[2], U[0,1], U[1,1], U[2,1], color='g', length=100, normalize=True, linewidth=2.5)\n", + "ax.quiver(points_mu[0], points_mu[1], points_mu[2], U[0,2], U[1,2], U[2,2], color='b', length=100, normalize=True, linewidth=2.5)\n", + "\n", + "# Display the rotated point cloud \n", + "ax = fig.add_subplot(1,2,2, projection='3d')\n", + "ax.scatter(points_rot[::subset,0], points_rot[::subset,1], points_rot[::subset,2], c='k', marker='.')\n", + "# Set the number of ticks on each axis\n", + "ax.xaxis.set_major_locator(MaxNLocator(nbins=5))\n", + "ax.yaxis.set_major_locator(MaxNLocator(nbins=5))\n", + "ax.zaxis.set_major_locator(MaxNLocator(nbins=2))\n", + "ax.set_xlabel(\"PC1\")\n", + "ax.set_ylabel(\"PC2\")\n", + "ax.set_zlabel(\"PC3\")\n", + "ax.set_aspect('equal')\n", + "ax.azim = 90\n", + "ax.elev = -45\n", + "ax.set_title('Rotated point cloud')\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Fit polynomial surface in PCA-space" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Choose 0 or 1 - make sure arrows point into cranium\n", + "normal_direction = 1\n", + "\n", + "popt, dz_dx_func, dz_dy_func = fit_polynomial_surface(points_rot, normal_direction=normal_direction)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Generate sample point coordinates in PCA-space" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "points_flat, grid_shape = generate_sample_points(points_rot, popt, dz_dx_func, dz_dy_func, normal_direction=normal_direction, extra_distance=15, min_steps=20, step_size=0.4)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Visualise the sample points and the original surface" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Transform into world coordinates\n", + "points_orig = points_flat @ U.T + points_mu\n", + "\n", + "fig = plt.figure()\n", + "ax = fig.add_subplot(111, projection='3d')\n", + "ax.scatter(points_orig[::1000, 0], points_orig[::1000, 1], points_orig[::1000, 2], alpha=0.3)\n", + "ax.scatter(points[::30,0], points[::30,1], points[::30,2])\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Sample points_orig's values in the original volume" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# The original .nii file of the cranium\n", + "path_nii = 'YOUR_FILE_PATH.nii'\n", + "volQ = sample_in_original_volume(points_orig, path_nii, grid_shape) # Can take some seconds" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Save the new sampled volume (Choose name and/or path)\n", + "name = 'CHOOSE_NAME.nii'\n", + "\n", + "niiPCA = nib.Nifti1Image(volQ, np.eye(4))\n", + "nib.save(niiPCA, name)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Suture extraction" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Load the NIfTI-file back in" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Load the .nii file in\n", + "nii_img = nib.load(name)\n", + "# Get the data as a numpy array\n", + "volQ = nii_img.get_fdata()\n", + "\n", + "# Print max and min values\n", + "print('Volume minimum: ', volQ.min())\n", + "print('Volume maximum: ', volQ.max())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Animate moving through the layers of the \"unfolded\" volume" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ani_through_volume(volQ)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Crop slices" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Use \"z\" to find slice of interest. Choose the layer with a clearly connected suture that is closest to the outside of the cranium.\n", + "\n", + "Crop using tool so that the suture does not connect onto outside air." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "interactive_plot_slice = interactive_plot_slicing(volQ)\n", + "interactive_plot_slice" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Remove unwanted sutures interactively" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib qt5\n", + "\n", + "z_vals = [60,61] # Choose z_vals based on the layer chosen before.\n", + "\n", + "new_volume = interactive_process_volume(volQ, interactive_plot_slice, z_vals)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Choose layers and thresholding value" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Choose thresholding value which connects the suture." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "interactive_plot_thresh = interactive_plot_threshold(new_volume)\n", + "interactive_plot_thresh" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib qt5\n", + "ani, points, suture, coords_orig = extract_suture(new_volume, points_orig, interactive_plot_slice, interactive_plot_thresh, z_vals, grid_shape)\n", + "HTML(ani.to_jshtml())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Skeletonize" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "- Chooses only ONE layer to skeletonize and represent the suture (curve)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Choose layer\n", + "z_val = 0\n", + "\n", + "# Create skeleton\n", + "skeleton = skeletonize(suture[:,:,z_val])\n", + "\n", + "# Get the 2D coordinates of the skeleton points\n", + "y, x = np.where(skeleton)\n", + "\n", + "# Create a 2D plot\n", + "fig, ax = plt.subplots()\n", + "ax.scatter(x, y, color='r')\n", + "\n", + "# Set labels and title\n", + "ax.set_xlabel('X')\n", + "ax.set_ylabel('Y')\n", + "ax.set_title('2D Skeleton')\n", + "\n", + "# Show the plot\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Get the 3D coordinates of the skeleton points\n", + "skeleton_indices = np.where(skeleton == 1)\n", + "skeleton_coords = coords_orig[skeleton_indices + (z_val,)]\n", + "skeleton_coords = np.squeeze(skeleton_coords)\n", + "\n", + "# Separate the x, y, and z coordinates\n", + "x, y, z = skeleton_coords[:, 0], skeleton_coords[:, 1], skeleton_coords[:, 2]\n", + "\n", + "# Create a 3D plot\n", + "fig = plt.figure()\n", + "ax = fig.add_subplot(111, projection='3d')\n", + "ax.scatter(x, y, z, color='r')\n", + "\n", + "# Set labels and title\n", + "ax.set_xlabel('X')\n", + "ax.set_ylabel('Y')\n", + "ax.set_zlabel('Z')\n", + "ax.set_title('3D Skeleton')\n", + "\n", + "# Show the plot\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Export the skeleton" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Choose name for skeleton\n", + "np.save(\"CHOOSE_NAME.npy\", skeleton_coords)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Curve Analysis" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The following functions will perform curve analysis as described in the method section of our rapport." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Get data from Blender (ply file)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "vertex_data = load_data_ply(\"OGAgamodon.ply\")\n", + "my_scatter(vertex_data)\n", + "len(vertex_data)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Calculate the mean points\n", + "* this depends on resolution chosen in Blender" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "mean_points = mean_points_func(vertex_data)\n", + "len(mean_points)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Sorting points\n", + "* Skip this step if the data came from Blender\n", + "* When points are extracted using the skeletonize method they have to be sorted, else the interpolation and smoothing function wont work.\n", + "* This function takes a bit to load (a few seconds)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sorted_points = interactive_sorting_points(skeleton_coords)\n", + "sorted_points" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Interpolation and subsamling of points\n", + "* The function uses the simple linear interpolation\n", + "* Note that if the data came from Blender \"sorted_points.result\" needs to changed to \"mean_points\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "interpolated_points = interactive_interp_points(sorted_points.result)\n", + "interpolated_points" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Smoothing\n", + "* Usually it works just fine with poly = 1, which means it fits a linear polynomie to the points.\n", + "* In generel the window length has to a relative higher value than the poly order, so keep that in mind.\n", + "* Also mind how much smoothing there is done (the window length) since the suture may loose its characteristics. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "smoothed_points = interactive_smooth_points(interpolated_points.result)\n", + "smoothed_points" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Calculate radii\n", + "* This is always done on the smoothed curve " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "centers, radii = calculate_radii(smoothed_points.result)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Finding the turning points\n", + "* This is just to get an idea of where the points are located this is NOT the actual amount of turning points" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "thresholded_points = interactive_thresholed_points(smoothed_points.result, radii)\n", + "thresholded_points" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Find the actual amount of turning points" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "turning_points = interactive_turning_points(smoothed_points.result, radii)\n", + "turning_points" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Finding the characteristics of the suture" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "curve_characteristics(smoothed_points.result)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Fagprojekt", + "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.10.14" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/Notebooks/Curve_Analysis_Script.ipynb b/Notebooks/Curve_Analysis_Script.ipynb new file mode 100644 index 0000000..e45843e --- /dev/null +++ b/Notebooks/Curve_Analysis_Script.ipynb @@ -0,0 +1,380 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Curve Analysis" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The following functions will perform curve analysis as described in the method section of our rapport. Note all the functions and packages are imported from \"myHelperFunctions\", which includes complete documentation of those said functions. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Loads the functions" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "starting functions\n", + "Functions loaded\n" + ] + } + ], + "source": [ + "from myHelperFunctions import*" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### get data from blender (ply file)" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaEAAAGTCAYAAACBGq1kAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAC7NUlEQVR4nOydd3gc9bX+3y3qvfferWo1WzK2gYAx1SYXSAi5DiQmQEJyg5PrAAaMTWgOPblAKKGEACGmhmaKG7bcZK1WXVr1XlcraXub3x/6zbC72r6zVfN5njwxq9nZ2d3Zeeec7znvYREEQYCBgYGBgcENsN19AAwMDAwMqxdGhBgYGBgY3AYjQgwMDAwMboMRIQYGBgYGt8GIEAMDAwOD22BEiIGBgYHBbTAixMDAwMDgNhgRYmBgYGBwG4wIMTAwMDC4DUaEGBgYGBjcBiNCDAwMDAxugxEhBgYGBga3wYgQAwMDA4PbYESIgYGBgcFtMCLEwMDAwOA2GBFiYGBgYHAbjAgxMDAwMLgNRoQYGBgYGNwGI0IMDAwMDG6DESEGBgYGBrfBiBADAwMDg9tgRIiBgYGBwW0wIsTAwMDA4DYYEWJgYGBgcBuMCDEwMDAwuA1GhBgYGBgY3AYjQgwMDAwMboMRIQYGBgYGt8GIEAMDAwOD22BEiIGBgYHBbTAixMDAwMDgNhgRYmBgYGBwG4wIMTAwMDC4DUaEGBgYGBjcBiNCDAwMDAxugxEhBgYGBga3wYgQg8shCAJqtRparRYEQbj7cBgYGNwI190HwLC60Gg0UKlUkEqlYLPZ4HA44HK54HK54HA4YLOZ+yIGhtUEi2BuRRlcABn9qNVqEAQBpVIJFotFRUMsFgssFgscDgd+fn7gcDiMKDEwrAIYEWJwOlqtFmq1GhqNhnpMqVTqCQxBECAIQi9Fx2az9aIkRpQYGHwPRoQYnAYpKiqVSi/aIR8DABaLZfK5jCgxMPg+jAgxOAXd9BsASoAAQKVSYW5uDmFhYfDz87N6f4aFDCwWS0+UuFyuSVFjYGDwTBgRYqAdMtLRaDRgsVh60YpYLEZzczNkMhm0Wi3CwsIQFRWFqKgoREREgMPhWPUalkSJFCZGlBgYPBtGhBhogyAIaDQaqvyazWZTIkAQBMbGxtDZ2Ym0tDSkp6dDpVJBJBJhfn4e8/PzUCqVCA8Pp0QpPDzcLlEiU39k9Z1uoQMjSgwMngUjQgy0QBAEFf0A+uk3tVqN9vZ2zM3NoaysDNHR0VR1nK5IyWQyPVFSq9WIiIhAZGQkJUrWrgGZEiXDNSVGlBgY3AsjQgwOo9VqoVQqV0Q/ALCwsAA+n4+goCCUlZUhICCA2l5XhAwhCAJSqRTz8/OUMGk0GkqQoqKiEBoaapUokac4I0oMDJ4HI0IMdkOm38jqN8P02+DgIHp7e5GTk4OsrCy9vykUCrMiZOy1JBIJFSWJRCIQBLFClKzZHyNKDAyeAyNCDHZhLv2mVCrR0tICsViM8vJyREVFrXiuUqmknmfv64vFYj1RYrFYeqIUEhJisyhptVrqOYwoMTA4H0aEGGzCVO8PydzcHFpaWhAZGYmSkhKjJdh0iJAhWq0WYrEYQqEQIpEIIpEIHA4HUVFRlDAFBwdbLUq6fUparRbDw8PIyspCQEAAVQ5umHpkYGCwHcY7jsFqzPX+aLVa9PX1YXBwEIWFhUhNTXXpBZrNZiM8PBzh4eHU8SwuLmJ+fh7T09Po7e0Fl8uloqTIyEgEBQUZPUbd98XhcKBUKjEyMoKMjAzI5XJqGzabDT8/Pz3fO0aUGBhsgxEhBqsgox+tVgsAegUBMpkMfD4farUa69evR1hYmLsOk4LNZiMyMhKRkZHIysqCRqOhRGliYgLd3d3w9/enRCkqKgqBgYFG96UrSBwORy9S0hUlY2asjCgxMJiHESEGs+im34xVv01NTaGtrQ0JCQkoKiqyuq/H1ZCpOXJ9SqPRYGFhAfPz8xgbG0NXVxcCAwP11pQCAgKM7sswUtJN3TGixMBgG4wIMZjEsPhA9yKq0WjQ3d2N8fFxFBcXIykpyZ2HajMcDgfR0dGIjo4GsNzLRK4ljYyMoKOjA8HBwYiKirIY2ZkTJYVCAblcDjabvaLQgRElBgamMIHBBLrWO4YXS7FYDD6fDzabjfLycgQHB9u0b2cUJtCNrpuDUCiEVCpFSEiI3pqSLb53pCiRGCsJZ0SJYTXCiBCDHpZ6f0jrnfT0dOTl5dnlYu0NIqSLUqnEiRMnsGbNGiqFJ5VKERoaqidKXK51iQVTomSYvrOlj4qBwVth0nEMFObSb7rWOxUVFYiLi3PnobqFuLg4JCYmAgAUCgUVKQkEAsjlcoSFhVFrSpGRkSbXx0hxIQWcFCW1Wg2VSqU34I8RJQZfhxEhBgDfj922ZL2zYcMGkwv2q4mAgAAkJCQgISEBACCXy6mm2e7ubigUCsqMNTIy0qxDuC2ixEydZfA1GBFa5ZAebSqVCgEBAVZb7zDoExgYiKSkJKpAQyaTUW4OExMTUKlUKxzCTYmItaLEDPhj8AUYEVrFkMUH/f39UCgUKCkp0bPeaW1txdLSEqqrq1dY7zgKi8WCLy9HBgUFISgoCMnJyZRDOClKo6Oj0Gg0iIiIoEQpLCzMZlFSqVTU2hozdZbBW2FEaBVirPeHtOAB9K13NmzYYHUVGINxWCwWgoODERwcjJSUFD2H8Pn5eQwPD+uZsUZGRiIsLMxk1GlJlJhIicGbYERolWGs+IDNZlMeaaT1TkFBAdLS0pj0mxNgsVgICQlBSEgIUlNTKTNWstBhcHAQAKw2YzUmSuRNhkqlwujoqF4FHzMKncGTYERoFWGq94fNZkOtVuPcuXNQqVQeY72zWmCxWAgLC0NYWBjS0tJAEASWlpYwPz+Pubk59Pf3UzZEpCiZM2MlixhIpqenERQUpBcpMVNnGTwFRoRWAebGbgPLzadzc3NISUnxaOud1QKLxaLMWDMyMqDVailRmpmZocxYdUXJlBkrsPz9k9EP+d/kDQkjSgzuhhEhH8da653Q0FCUlJS481AZTMBmsxEREYGIiAhkZmZCq9ViYWEBIpEIU1NT6OnpocxYSWEKCgqinq+73gfoR0q6s5QMRYmZpcTgChgR8mHM9f7oWu/k5uZibm7OjUfqHXjKRZjNZlMREOkQTjo5kA7hAQEB1Da6g/oM0fW8A/RFSalUUhNwGVFicBaMCPkgunN/rLHemZiY0LOQYfAujJmxkqI0MjICuVyOzs5OxMTEUMLk7+9vdF/mREmhUJgtCWdEicEeGBHyMbRaLdRqtU3WO2SJNoNvwOVyERMTg5iYGADAyZMnkZqaCqVSiaGhIbS3tyMkJERvTclUGb65WUqGokSuJzFTZxlsgREhH8HS2G1z1jvuahw1XKtgcB6kdRCg7xA+MDCAtrY2yoyVHARoTpSMja0wHPDHTJ1lsBZGhHwAc2O3rbHeIfuEGHwTQ7H38/NDXFwcFQkrlUrK966vrw9SqRRhYWFUlBQREWHSIZwRJQZHYUTIy9Ht/dFtWASst97xdQud1Y6liNPf31/PjFWhUFBuDqQZq6EoWTJjBcyLEjN1loGEESEvxVLvjy3WOywWi4mEfBhbbzACAgKQmJhIja2QyWRU+q6zsxNKpRIRERHUmlJERIRF3zvA9Ch0Zurs6oYRIS/EXO+PPdY7TGGCb+Po2htpxpqUlKRnxioSiTA+Pg61Wm2zGSugL0oajQYajYYZhb4KYUTIyyD7N4xFPzKZDC0tLVAqlTZZ7zDpON+GzgIQS2asIyMj0Gq1eqIUGhpqs0M4KUpkBZ5arUZ0dLSe7x0jSr4BI0Jegrmx2wAwNTWFtrY2JCQkoLq62ibrHaYwwbdxZhWiMTNWiURCidLQ0JCeQzgpSraYsQqFQszMzFDWRMzUWd+CESEvwFrrneLiYmqomi0wkZB5vP2zcWUpPIvFQmhoKEJDQykzVrFYTInSwMAAWCyWnsWQJYdw0vvOz8+PmTrrgzAi5OHoRj+Gd3u61jv19fUIDg626zVcXZjA3LW6DlJA3fV56zqEp6en65mxzs3Noa+vDxwOh4qSIiMjVziEk6lncn/M1FnfghEhD8VS74+h9Y4jPzKmMMF3cbcIGaJrxgosC8zi4iLm5+cxNTUFgUAALpdLiVJUVBQ1esQYzNRZ74cRIQ9Ed+opAL0fjCnrHUdg0nG+j6eIkCHknKTIyEjKjJUUJdKMlc1mw9/fHxMTE4iKikJgYKDJ/TFTZ70PRoQ8CGNjt01Z79TX15v9MdoCKUKMjY7v4WmRkCV0U3PAcjq6o6MDcrkcY2Nj6OrqQmBgoF6hg64FlSGWps6S2+iKEjN11rUwIuQhmCs+IAgCQ0NDEAgEyM7ORnZ2Nq0/Et0fqCt+fFqtFsPDwwCAmJgYswPZGBzD20TIEA6Hg4CAAAQGBiIvLw9qtRoikQgikQgjIyPo6OhAcHCw3pqSKYdwYOXUWXMD/gyr7xicAyNCHoCpsduA9dY7jqDb6OrstIRcLkdzczNUKhX8/PzQ29tLDWSLjo62eGfLYBveLkLA8nlJOn5wuVzExsYiNjYWgL4Z6+DgIMRiMUJCQvREyZJbiLWixEyddQ6MCLkROq13HEE3EnIms7Oz4PP5iI+PR15eHvWaCwsLEAqF1J0teRGJjo5GZGSkSfNMBsv4wlqfuZsjY2aspCjpmrGS6TtL55O1osTMUqIP5tftJui23nEE3bSfMyAIAr29vRgcHERRURFSU1Op9244kE2lUlE9JQKBAHK5HOHh4XrmmcwisvX4QiRENmdbg7+/P+Lj4xEfHw9g2YyVFCXyfCLNWMmCCHON3bqixEyddQ6MCLkBc2O37bXecQRd8aMbhUKBlpYWyGQyq96Pn5+f3kVELpdDKBRifn4eY2Nj0Gg0iIyMpFJ35rrvGXxDhBxJEwcEBOg5hMvlcsr3jnQI173JCQ8PN+sQDjCj0OmGESEXYtj7Q6f1jiM4KxKan59Hc3MzoqKisHbtWrvSaoGBgUhOTkZycrKeJYxQKMTAwADYbDZ1AYmOjkZQUBCt78Hb8RURouv4AwMDkZSURDmLkGasZEm4SqVaIUrmepQAfVEyNXWWESXTMCLkIrRaLeWBlZmZqdd8Sof1jiOQx0JXJKQ7SC8/Px/p6em0/OgMLWF0Gx0nJyfR09ODgIAAKkqKiooyWyllz+t7G75Qdu/MghnSIZy8ydEVpdHRUWg0GocdwklRam9vR0JCAlUswTiEL8OIkJPRXdiUSCQYHx9HdnY29Xe6rHcchS7XBJVKhdbWViwuLqKmpgaRkZGOH5wJDBsdyfJd0jizvb0doaGhlChZyv/7IowIWY8lh/Dh4WE9M9bIyEiEhYVZNGMFlkVpYWEBcXFxzNRZAxgRciKGxQdcLpeKNgiCwPj4ODo6Omix3nEUOlwTFhYW0NzcjNDQUNTX19MahViDYfkuObZaKBRS+X/yrjY6OtrsXa2v4AsiZEthAp0YcwgXi8V6JeEA9BpnzZmxajQavZQcM3V2GUaEnISx3h8OhwOtVgu1Wo2Ojg7Mzs7SZr3jKI6k4wiCwMjICLq7u5GTk4OsrCyLPxpX/Kh0x1YbplpGRkZAEITeepKhcaYv4Asi5KpIyBK6ZqykQ7iuGWt/fz8VnZPnle45RVaD6u7PWPrO2NTZjz76CBUVFSgrK3P9G3cyjAjRjLneHzabDbVajYaGBtqtdxzF3nQc6WUnFApRVVVFlVp7GsZSLWKxGEKhkHJzJo0zyfSdp3w3juILIuSJ74HFYiE8PBzh4eHIyMjQcwifmZlBb28vuFwuJUrmjFjJ/ZkSpb/+9a/47W9/y4gQg3ksWe9MTk5CrVYjKyuLdusdR7EnHbe0tITm5mYEBASgvr7eq5wOdO9qyQsI2TRLepQFBQVRi9HeChMJuQ5dh/DMzEzqnBKJRJiamgJBEODxeFQTdlRUlNlqTl1RkkgkCA0NddVbcSmMCNGEud4f0npncXERAJCTk+OuwzSJrdNVx8bG0NHRgczMTOTm5tp8ofO0C6NuqTegbwczNDQEAGhsbKSipIiICK8ocmBEyH3onlNpaWk4fvw48vPzsbS0RDmEBwQE6I2tMHYjRxZIePPNkDm875v1MHRt4k1Z75w8eRJsNhvV1dXUczwNayMhjUaDtrY2dHV1oaKiAnl5eV5/kTMGaQeTn5+PyspKAEBqaioUCgU6Ojrw3XffgcfjYWhoCIuLix75nQK+IULuKkygEzI7Ehsbi5ycHFRVVWHjxo0oKCiAn58fRkZGcPLkSZw+fRpdXV2Ympqi+oyA5UgoJCTEptf88MMPcemll1ImwVlZWbjxxhsxMjKit93i4iJ27dqFjIwMBAQEICMjA7t27aJumo3x9ttvo7a2lrLYuuKKK9DY2GjT8ZEwkZADkEUG1lrvkNbxWq3W4+6irSlMkEgkaG5upsrJV1tjaGJiIpKSkvRKd4VCIRUpGTbNesLF31PF0Ra8NRLShVzX0n0fXC4XMTExiImJAQCjLQbvvfce1Go1FhYWqCZ3SxAEgdtvvx0vvfQScnJy8OMf/xhhYWEYHx/HsWPHMDQ0hLS0NADLv+nNmzejubkZl156KW688Ubw+Xw8/fTTOHLkCE6cOLFC/B555BHs2bMH6enpuP322yEWi/Huu+9iw4YNOHToEC688EKbPhtGhOxAt/fH2NhtU9Y75AnoiSJkqTBhcnISbW1tSElJQUFBgcMXBXLR1RsxLN3VXZCenp6GQCCAv78/lbqLjo52ebk6iS9EEZ5amGALhpVxxjDlEP75559DqVRi8+bNKC0txUUXXYTrr78e9fX1Rvfzl7/8BS+99BJ+/etf49lnn13xurpiduDAATQ3N2P37t14/PHHqcf37t2L/fv348CBA9i3bx/1uEAgwN69e5Gfn4+zZ89SE3J/+9vfora2Fjt37kRXV5dN7ijefXa6AdJ6R6lUGhWgqakpNDQ0ICQkBHV1dXp5XF0R8jRMpeO0Wi06OzvR1taGkpISFBUVef1FjW7IBenMzExUVlZi06ZNKCoqotIsJ06cwJkzZyAQCDA7O2v1HS0deKvQk5A3K95+zlmqjDOGn58ffvjDH+Kpp54CALS1teHee++FTCZDa2ur0efIZDLs27cP2dnZeOaZZ4wKHykQBEHglVdeQWhoKB544AG9be655x5ERUXh1Vdf1TuHXnvtNajVauzZs4cSIAAoLi7Gjh070NfXh8OHD9v0PplIyAZ0e38MQ2trrHfIdB2ZvvMkjBUmyGQyNDc3Q6vVoq6uzuactK9g64XclDO4UChc4QweHR1t1p/MUbx9TcjYiHuS4TkJ/nluDAF+LNy5McNt0aY1WBMJmUIikQAAMjMzUVRUhBtuuMHktl9//TWEQiFuvvlmaDQafPLJJ+jp6UFkZCQuueQS5ObmUtsKBAKMj4/jsssuW/HbDgwMxKZNm/Dxxx+jt7eXGr1y9OhRAMCWLVtWvPZll12GF198EceOHTP6d1MwImQFlub+2GK9Y2sVmqswjIRmZmbQ0tKChIQEFBUVeVz60JswdAbXbZptbW2FVqulSnajo6PNdt3bii+KkEqlwg9f4WFQKKcee+30BACgOikAr/18nWsP0gocScFLJBL4+/tb1QJBFgdwuVyUl5eju7ub+hubzcZdd92FJ554AsCyCAGgBMYQ8nGBQKD379DQUCQmJprd3hYYEbKApd4fW613PFmEtFotCIKAQCDA0NAQ1qxZg5SUFHcfms9haJopkUiocRW6zuDkmpIjBSDeLkLkjZHu7+qX77TrCZAujRMKlD58HEnBwFd3bXLJMVqDPek4EnJarDXf4/T0NADgySefRGVlJc6ePYuioiLweDz88pe/xJNPPomcnBzccccdWFhYAAC9tJou4eHhAEBtR/6bvJmyZntrYETIDOScEGPRj73WO54qQmw2GyqVCufOnYNCoXDZLKPVjq4zeHp6OuUMLhQKqV6SwMBAPVGyZcKut4sQ+VvRfQ+NI6ZLh0kmpEDpw8dx29pg3L6lwu3TeR1Nx1lrbEx+Xv7+/vjoo4+QnJwMANi4cSMOHjyIsrIyPPnkk7jjjjvsOhZnwIiQEcj0G1n9ZihACwsL4PP5dlnveKoIqVQq9PX1IT4+HpWVlW7/0XoirriY6zqDA/pluwMDA2hra6Mmg1rjDO4LIqT7+1PZuJ76N54Uf+M16D2WE83Bh7fXu/RzcVSErB3eSEY11dXVlACRFBcXIzs7G729vRCJRNS2piIXsk9IN1KKiIiwaXtrYK40BlhKvw0NDUEgECA7O9su6x02m+1RhQkEQaC/v58Ks8vKylzy4/TmC6MrMeYMTqbuurq6oFQqERERQUVJhs7gviBCusfvR8PaZJ9Qg7JHvkN1aihe+1mlw/uzBkfXhKwtCiooKAAAkyNUyMdlMpnFNRxja0Z5eXk4deoUJicnV6wLWVpjMgUjQjroRj+Gpdek9c7S0hKqq6spexdbIZ20PQHyPYnFYsTExCAiIsKrL1irAX9/fyQmJiIxMVHPGVwoFFLzbnSbZr29x8ZYo2paVCBG5o2vCdlC46gYpQ8fx5v/XYK16c413nVkTcgW37iLLroIANDZ2bnibyqVCr29vQgJCUFcXBwSExORnJyMkydPrhA6uVyO48ePIzk5Wa+ibvPmzTh16hS++uor7NixQ2//hw4doraxBe8uvqcJXesdYwKka72zYcMGuwUI8Jx0nEgkQkNDA1gsFjX7xxOOi8F6SGfwlJQUlJaWYuPGjVi7di0iIiIwOzuLc+fOoa+vDxKJBBMTE1AoFO4+ZJsx1iP04a1rab1w7fhHGy77y0ka97gSV60J5eTkYMuWLejt7cUrr7yi97fHHnsMIpEI1157LbhcLlgsFnbu3AmxWIz9+/frbfvoo49ifn4eO3fu1LsW3nLLLeByuXj44Yf10nLt7e148803kZOTg4svvtim97fqIyGy98dYKagx6x1H7yrdLUIEQWB4eBg9PT3Izc2lRo3TNVnVWrz57txTMRwtoNFo0N/fj9nZWYyNjaGzsxPBwcF64889fe3PWCQU4OcH/p5NuPuDVnzWOU/L64wvalD68HF8+atKpETR71at1Wrt/qxtddB+/vnnUV9fj1tvvRUfffQRCgsLwePxcPjwYWRkZODPf/4zte3u3bvxySef4MCBA+DxeKiqqgKfz8cXX3yBiooK7N69W2/f+fn5ePDBB3HfffehrKwM1113HSQSCd555x2oVCq8/PLLNr9Pzz4DnYiu9Y6x6jdT1juO4k4RUqvVaGtrw/z8/IqUoiND7Rg8Ew6Hg6CgIISEhKCsrIyygREKhejr64NMJqOKHKKjoxEREeFxzgTmfOMe+2EpHsNyCXPds020vN7W55uwrSQWf9q2hpb9kWg0GrtHnYjFYptEKCcnB42NjXjggQfw5Zdf4quvvkJiYiJ+/etf44EHHtArsQ4JCcHRo0exb98+HDx4EEePHkViYiLuuusu7N271+ha1J49e5CZmYlnnnkGL7zwAvz9/VFfX4/9+/ejpqbG5vfHIrzd18MODIsPDNNv09PTaG1tRUJCAgoLC2m9W2xqakJ0dDQyMzNp26c1LC0tgcfjISgoCGVlZSt+EJ2dnWCxWCgsLHTJ8ZDl75520TOGXC5HQ0ODzWkGT2BkZAQikQilpaUr/iaXy6mmWaFQCLVaTRU5REdHW12R5UxmZmYwMDCA2tpai9ve+W4zjvVZLt+2ltY99PUZtba2IiIiAunp6TY/9+677wYA/N///R9tx+NJrLpIyNjYbRJrrHccxdVpLwAYHR1FZ2cnsrKykJOTY/TCYs9QOwbPx9x3GhgYiKSkJD1ncLLybnBwkCoX122adbUo2eKg/dcfVwAAvmkbx72f9ULmoEVf6cPH8d2u9YgMctwOyNE1oYSEBIePwVNZNSJEp/WOI7iyRFuj0aCjowPT09NYu3YtVeZrDE/1tGNwDGtLtHWdwdPS0ihncKFQiKmpKfT09OgNYHOVM7g9YxwuKUnGJSX6PTJarRa1j52Awsb7rI1PnQYAbMgKx1+uL7apUVgXOvqEfJVVIULmen+A76eEWmu94wiuKtEWi8Vobm4Gl8vFhg0bLDbUuiNCY3A+9vYJ6Y6qzsrKgkajoZpmh4eH0dHRgdDQUL2mWWcUOdDloM1ms9F47yZ0T8zjur8bd6A2x8mBRVQeOIXEEDbev70W4YG2CbAjfUJSqZQRIW/G3Nhte613HMEVhQkTExNoa2tDWloa8vPzrfoRM4UJvgldzaocDkdvAJtSqaSKHHp6eqBQKBAeHk6l7uhyBqd7oF1BUhRa92zC5X9pwOii7fm6SYkWG548jeggNj69rRJhIdZlTOjwjvNVfFaEyLk/5OwWOq13HMGZ6TitVouuri6Mj4+jrKzMpjwysybkmzjLMcHf39+oM/jM7CxGRkZAEAQtzuDOarb94jf1ONQ+jj981GvX84UyLeqfaURKuD8O7lyL0CDzlW+OpuN82cfRJ0XIsPdHt/qNDusdRyCNQulGKpWCz+eDIAi71rSYdJxpvPlzcYVtD0EQqH3qnMGjLDy4SYFczKO/vx9cLldvPcnamz5njva+rDgZF+bGoPqJM3bvY2xRibqnziA1lI2Pf7UO/ibWjFxl2+ON+JQIWRq7TZf1jiM4IxIiS8oTExNRWFho18nOpON8E1eIUNkj3xl9/MHjcgByLBuzaAHM/f//EQD0j+nlHxVife7KEQHOFCEACAgIQOueTSh9+LhD+xkVa1F14BQSglj4ZtfGFX+3NxIiqxaZNSEvgCAIiMViTE5OIjU1dYUACYVC8Pl8REZGYsOGDXZXuTgKnYUJWq0WAoEAw8PDKC4uXuGaawtMJOS7OFOE/nHCtgFmyxgeD4Fb/9UFoEvv0etKo/HT4iCX9JLRIUQAMCUjUPrwcbTcu1Ev++KImPp6dZzndwpaAdn4KJVK0dPTo7f+Q16oz58/j5ycHFRUVLhNgAD6ChPkcjnOnTuHmZkZ1NXVOSRAABMJ+SrOjoQOHJugYS8rRQkADrYKsf3dMfz44BSUSiUNr2Oe1j2bwL/nAlySF+nwvip1okMy88GUaBvHqyMhw94fLperl+pylvWOI9AhQnNzc+Dz+YiNjUVVVRUtpbGuLkyYm5vD1NSU13iYeSt0lTi7Fl1RIqAGUPXn0+AAOHpnOSLCw50mrGw2G0/fUEb999unh/Hot4M270e37o78vdsjQhqNBnK5nBEhT8RY7w+Xy6VC39nZWcp6h64LNR04IkLk7J/+/n4UFhZSaUe6jssVIqT7HuLi4igPM7K8Nzo6GuFOvMisNpwZCZ0bEjplv/roOJoA2PjXZoSxWXhlWyxV5GCvJ5s1/GR9On6yPh1ft41j18e2VdK9dEyAX27Og0ajWbE8YC1isRgAGBHyNEz1/pB3Gp2dnZiYmHCa9Y4j2CtCSqUSLS0tkEqlWLduHTXPnS5ckY5TqVRUYUhNTQ0CAwPBZrP1ZuKMjIwAAHWBiY6ORlBQkFOPy5dx5o3Fz99qc9q+TcPCkhb40Yez2Jo2i8tT1QgJCaEiamdF1ZeWJKO1JBn/Pj+K/V/2W/Wcj1qnKREybBGxFolEAoARIY9Bt/fH2NhtuXx50JVIJHKq9Y4j2FMdNz8/Dz6fj4iICNTV1TllTcvZkZBYLEZTUxOCg4NRV1cHLpdL5fmDgoIQFBSE5ORkEASBpaUlzM3NYXJyEj09PQgKCqIEyVmd+b6Kt09WNceXI8CXI1z857+ToJDJ9JzByaZZup3Br69KxfVVqXji6268cXbK7LZsYvl37miPUGBgoE+f817zzrRaLdRqtUXrHQAoLy/3SAECbKuO0+1pysvLQ0ZGhtMuKM6MhCYnJ9Ha2oqMjAzk5eWZfS3dmThZWVlQq9VUlCQQCCCXy/WcnsPCwnz2IksHhiL0r7PD+NPXg3rbPHxVNq4pT7Vpv5c+c4KOw6OFq/8xCA6Ac3+sh0bnfBkfH4dardZrmqXLGTw50LJb9y82ZACgp0fIl89xjxchS70/htY7ZMOmp2JtOk6lUqGtrQ0LCwuoqakxOTOeLpxRmEBWJo6MjNjs4EDC5XIRFxdHWSrJZDIIhUJqnDUASpBsaYK0FW+9COiK0C/fPo9TA5IV2+z5tB97Pu0HG8CRX1cgOtJyqndSYvkcvmFtHG6oiEd+UrTRz+8/TYO494thy2/CCjQAKh9vAAAc/HkZiouLQRAEJBIJNa6CdAbXbZq1N9X76HGZxW22r10e2+BIJCQWiz32hpouPFqEDK13DAXImPUOh8OhtvdErBGhxcVF8Hg8hISEUKO3XXFcdIqQUqlEc3MzFAoF1q9fT1tOOygoCCkpKUhJSdFzep6YmEB3dzeTujOAFKHPWieMCpAuWgCb/68ZAMC/5wKTaazrXmyw+LppYRzcf0WR2W2urszE1ZWZeo8pFAr84K+NWJDb39B93d9bAAB3/yAdN63PRGhoqEVncDJ9R+dvjbxWOeIbR5Zne+tNkDV47C9Ud+4POX6axJz1DofD8eiRBOZEiCAIjI6Ooqury+WWQnSm4xYWFsDj8RAREYHKykqnCYGh07Pu5FDSVHO1p+5IEbr7E9uaSssfPWFyqFv3nOWbvM9/u8Gm1yMJCAjAid8vP/eGl8+hc9pyxGGKx74dxmPfDuOO+mT86qJcs87gQ0NDaG9vp5zByUmzxs7db5qtK0wgcXRNyJctewAPFCFLc38sWe+4alSCvZgqTNBNK1ZWVlJuxa48LjoiIXKAXm5uLjIzM1160ffz89NL3ZFD2oRCIYaGhqhUTHR0NGJiYpxa2uspEASB33w8ZNdz/9M0hKsrM/Qe651csPg8ur7x926twdmzZ3HXETnsMLymeKFhHC80jOPo/1QjJvT71JYxZ3Ayddfd3a3nDE7exLDZbNz12ajF16xM+F50GN8483iUCFma+2ON9Y6np+M4HA4IgtDL1ZOzf/z8/Fzq6K2Lo2tCWq0WHR0dmJqacouIGiM4OBjBwcFITU2FVqvF4uIitWDd1dWFkJAQvdSdvRcKT2ZeqkL/on3f671frBSh61/lW3ze2d9bHsVtLQRB4JOdxeAGBeOCpw1NUm3jwmcbAQC7NifhlgvyVvzd398fCQkJ1Noluf44Pz+P0dFRaLVa/OGkdRL7xs7vI0FH03GMCLkI0nrHWPSj1WrR19eHwcFBFBQUIC0tzeQdtjek44Dv747Gx8fR3t6OjIwM5Obmuq273ZEmWrlcDh6PRzl4e2JfDzmqOjIyEtnZ2VCpVFQVFXnXS46yprOKyt3sazC/DmSJbc834ONf1VP/bc3t3Redc7h2bYpDr0tCXg8igoPQumcT7vu4HR+3zTm0z6eOTeCpYxMm040kuuuPGo0G6584CTml5ytNWEk2G7x1ZqqqedwuQmT6jax+M9b7w+fzrbbe8RYRUqlU6OrqwuTkJMrLy6m5LO7C3kiItBCKi4vDmjVrbPqxufMi7+fnR83DIQhCr+qOrKIiBclbq5NmxArM2L+kAgDon1dDIlMgxMK8HF0e+LwP28qTnDLU7k/bivGnbcAPnjmBaSsq9MxBGpZaEiORTIWNT50yeNSY393yYw9cUQiVSkVlahgRMo9bRchS+o0cUWCL9Y63iFBjYyM4HI7HRA62FiboFofYayFEip67Iw4Wi2U0dTc3N4fR0VEsLS0BAAQCgVel7m5+k0fLftY/dcbihdoQc4UNtmBqqN23v7sAAPCnz7vwL960Q69R+vBx/LYuFrdevGbF394/P4IHvxywYi/Lx7ghDhgYGEBbWxvVNCuVSu2+kZFIJF57E2QtbhMh3ejHsPRaq9Wiu7sbY2NjNlvveLoITU8v/2AiIiJQXFzsMeaSthQmqNVqtLW1YX5+3iU9TK5GN3WXk5ODxcVFNDY2Qq1Wo6urCyqVChEREYiJiXFoaqgzUalUGJ6nz3m64uHjaN6zCff/IA0PfTti1XNKHz6Ok7+vQ3ig/Q4flkYg3HdFIe67ohD8ESF++qb9NkLPnZrFc6eO47NfliE9LhIKpRKbnz4NiQ3Ly34s4MVfLguvQqGg0r3k/4vFYqowxtpKTYlEgpQUelKbnopbRIjs/zEmQBKJBHw+HywWyy7rHU8VIa1Wi56eHoyMjIDFYiErK8tjBAj4PhqxZPMikUjA4/GoIorVUGHm5+cHFouFoqIiasgYmbojp4bqNsy6oq/LEr//qMviNikhwJiVS0YaABc99R2O7NpotQgBwIYnTzkUEVnrAl6eFo3WPZswuyTFRc812v16V77UYtfzgrjA2T9+/z4DAgKQmJiIxMREqNVqBAcHIzAwEPPz81STtWHTrLHf3WpIx7nlKkj2/Riz3mloaEB0dDTWrVtnVxjqiSIkl8tx9uxZzM7Oor6+Hlwu1+PKyHULJkwxPT2NU6dOITY2FjU1NatCgAxhsVgICQlBWloaysvLsWnTJhQXFyMgIAAjIyM4ceIEzp49i97eXgiFQredi0d65s3+3Y8NfPm7TQixIYCblRG46dUzOPP7dTYdy7tnB23ansSeYXCxYcG0pAFtITXST0+ADCEIAkFBQUhNTUVpaSk2btyIiooKhIWFYWZmBmfPnkVDQwM6OjowOTkJhUJBPdeWqapkS4Sx/91+++0rtl9cXMSuXbuQkZGBgIAAZGRkYNeuXVhcNG1J9Pbbb6O2tpYyjb3iiivQ2Gi/6ANuTMeZs94h+zzswdNEaHZ2Fnw+H/Hx8dTCPV2D7ehENxIyhCAI9Pb2YnBwECUlJR7nTO5OdG1gcnJyqF4ToVCIzs5OqFQqvao7V6TuJhYsVyP86eoCAMDpe5cvnvUPH8eSFftumVTg9bNjNk0iffjrYfy4NtOqbXUhz0V7Mgb/vqUI17/WafPzbKUiOQT/uKXK7DaGhQm6/oiZmZnQaDRYWFjA/Pw8RkZG0NHRAR6Ph+7ubgwMDNj0/iMiIvC73/1uxePV1dV6/y2RSLB582Y0Nzfj0ksvxY033gg+n4+nn34aR44cwYkTJ1aUhj/yyCPYs2cP0tPTcfvtt0MsFuPdd9/Fhg0bcOjQIVx44YVWH6cubq+OI613AgMDaemR8ZRmVd0Ld1FREVJTvzeI9GQRMjwulUqFlpYWSCQSWgcDetoaijlsqRrU7TVxV+rud/9qNft3FoArSvR9/Br2bMLVL5zFoFBucf8vfDeK7WVJNglR6cPHbY5QyHPRnnOlMDkOrXvi8LcjPXihYRLOuC29oTIe919eaHE7S31CHA6HOh9ycnKgUqmoddeBgQH87ne/wzvvvINLLrkEl19+Oerq6kzuKzIyEg8++KDFYzpw4ACam5uxe/duPP7449Tje/fuxf79+3HgwAHs27ePelwgEGDv3r3Iz8/H2bNnERERAQD47W9/i9raWuzcuRNdXV12uaO4dVFicHAQZ8+eRUpKCjVfxlE8oVlVoVCgsbERExMTWL9+vZ4AAZ4pQuSPRPeCu7S0hIaGZZ+wuro62ifTepMQ2YOx1N2aNWvg7++P4eFhnDhxAufOnUNfXx/m5+dpOyc6ZswLyZUlxjMN/7mjFgEc676Ty/5vuXF0fbr16xWlDx+HTCfVZAny83Bk7fS2i/LRvGcTmu/eAD8ar3Z/u77IKgECbC/R9vPzw5YtW/Daa68hKysLL7/8Mm699Vb09/fjtddes/eQKQiCwCuvvILQ0FA88MADen+75557EBUVhVdffVXvWvDaa69BrVZjz549lAABQHFxMXbs2IG+vj4cPnzYruNxWyTU39+PwcFBo9Y7juDudNz8/Dyam5sRFRWFtWvXGr0z8EQRMkzHkU20WVlZyMnJ8XnBcAW6vUfAsk0MGSW1t7dDrVbrDfMLDg62+XMfFlquNNh/Ra7JvzXevdHq6Oa2t5vx8n9X4pKnT2JKat1vrvaJM3jzp8VYm2HZUcORdJwhHA4HTfcsR2InOsZxx4e2TUklCfIDGnbV23THb2+fEOkCnp2djQsvvBC33HKLxecoFAq88cYbGBsbQ1RUFOrr61FeXq63jUAgwPj4OC677LIVKbfAwEBs2rQJH3/8MXp7e5GXt+wscfToUQDAli1bVrzmZZddhhdffBHHjh0z+ndLuE2E0tPTkZSURHs6wl3pOIIgMDg4iN7eXuTn5yM9Pd3kBcRTRYjFYkGtVqOzsxNjY2Me0UTry/j7+1MVVOQFRygUYm5uDn19ffDz86MEyVqH5ye/Nm9UymbB4lDE5rs3oOKxkxZfq2FgeQH7m7s2YOtfTmFsUWXxOQCw4612HLg6B5eXmS89JnuE6L4BumBNMu6YXsQLJ23rL/r9ham4eUO2za/nSu+4yclJ3HzzzXqPbd26Ff/4xz8QGxsLYFmEAFACYwj5ODnHjPx3aGgoEhMTzW5vD25Lx/n7+zslH+6OdJxKpQKPx8PQ0BBqamosDp9zd7RmjpaWFszNzaG+vp4RIBfCYrEQGhqK9PR0VFRUYOPGjSgqKoKfnx+GhoasTt0d7jU/bK082XL6jMPh4PTvTa876PKzvy+n5b78jXXbk+z+Tx/G5s2XQthaGWcLv7qwEN/8qtzyhgBuXpeC1j2b7BIgwHHvOGvT4D//+c9x9OhRzMzMYHFxEadPn8bll1+OL7/8Etdccw0VWS4sLJvQ6qbVdAkPD9fbjvy3LdvbgtsLE+jG1Rf4hYUFNDc3IzQ01OrZP54YCc3Pz4MgCPj7+6OiomLVz+FxN7qL1cBymkU3dafRaPSq7qxN3T1+rfkZPyQhgX54/oYi/Oo98xVmTRMy1D9yHA33bgL/ngtQ/qj1E1e3Ps8zW6xgyi2BLhKiIvReX6FU4sPmSQwIFajPCsem/HiHX58sM7cnElKpVFAqlVaLkOH6zrp16/Dpp59i8+bNOHHiBD7//HNceeWVNh+Hs/G5K42r0nEEQWBkZATd3d3IyclBVlaW1SesJ4mQ7vtgs9nIy8tjBMgDCQgIQFJSEpKSkkAQBMRiMYRCIWZnZ/VSd5ZIirDeImpjXhwiAnqwoDB/U7dEAHe8cQ4v/KwGLfduRNkj31n9GlK5HMEmCpKcGQkZI8DfHz+uTad1n+Tv3B4REovFAOCQizabzcYtt9yCEydO4OTJk7jyyiupiMZU5EL2CelGPhERETZtb9Mx2vUsD8YV6Ti1Wg0+n4++vj5UVVXZPHzOU0RIo9Ggra0Nvb29VnvzMbgfFouFsLAwZGRkYO3atdi4cSMKCwt1vj/6JuSe+IN1w+lOjMqoyKV1zyaEWOnUs+7Jsyh9+DgqjBRDuFqEnAGZlbFHhCSS5SITR0c5kGtBUqkUgOU1HGNrRnl5eRCLxZicnLRqe1tw2zfsrDCbnNfjrIv80tISTp06BaVSifr6eqvuPg3xBBGSyWQ4c+YMxGIx9T7onK7qa3hydSA5nC03l6x60z1WArqiJJVKaR3jrsvmp79PxZ3evQnVKdZHXRosl3Df80E79Zi1lj2eDClC9pw/pHmpo2a5Z86cAbDsqAAsi0VycjJOnjxJCR2JXC7H8ePHkZycrHM+AZs3bwYAfPXVVyv2f+jQIb1tbMW7v2EjkF+YMy6mY2NjOH36NBITEx2yrTE1XdVVzM7OoqGhAREREVi3bh3Vn0XXdFUG92D6Qvf942fOnMGpU6fQ1dWF6elpqFSWK9pyY6wrIBLJgeGZ74sNXru5xqrn6fJp5xxKHz4OrVbrM5EQh8NxSISseW5HRwdEItGKx0+cOIGnnnoKAQEB+OEPfwhg+TzZuXMnxGIx9u/fr7f9o48+ivn5eezcuVPvdW+55RZwuVw8/PDDemm59vZ2vPnmm8jJycHFF19s83sEfHRNCFhOmdGVXtJoNOjs7MTU1JTDtkKAe8vIBwYG0NfXt8LFAbB9nAOD9/FWnz8OXFsAoVCIgYEBtLe3UyMHYmJiqBHWunx4+3qre4eufMl8sYG1lD96As9enQBCKsf5gRmUpUZaLC33RFxVnv3ee+/hwIED+MEPfoDMzEwEBASgra0NX331FdhsNl588UWkp3+/3rV792588sknOHDgAHg8HqqqqsDn8/HFF1+goqICu3fv1tt/fn4+HnzwQdx3330oKyvDddddB4lEgnfeeQcqlQovv/yy3ddbj/COoxPSFJWuSEMikaC5uZnW2T/uiITUajVaW1uxsLCA2tpao4uITCTkPYzOL6JBMI+CxFCUpUVb/Xs6O6VBTEwMNX5dLpdTVXdjY2PQarV6tkLk+X7+f9ej6s+nrXqNvx4R4M6L7Fsf0OV//jO1/I+T31fofXVbGZJiIx3et6twpDxbLBZb7TV40UUXobOzE01NTTh27BjkcjkSEhLwox/9CHfddRdqa/VHroeEhODo0aPYt28fDh48iKNHjyIxMRF33XUX9u7da1T89uzZg8zMTDzzzDN44YUX4O/vj/r6euzfvx81NbZHvSRujYTsneZpCbrKtCcnJ9HW1oaUlBQUFBTQlhpgs9lWpUHoQiwWg8fjUf58psrInfV9mMLVr+ftbH32O4yJTX9eX/+qwqr9iMQSRIYuX2QCAwORnJyM5ORkEASBpaUlCIVCTE1NoaenB4GBgZQg8f5Yj+1/a8SQyPycor81TFAidOfGZPz1u3Hr3qAVbPlbC879vpYWiy9X4MhUVVsctDdv3mzzmkxERASeeuopPPXUU1Y/56abbsJNN91k0+tYwrsTriZwVIS0Wi06OzvR1taGkpISFBUV0ZqbdmVhwtTUFE6fPo34+HhUV1eb7WNi0nGeyc2vNqD04eNmBQgALn2+GRWJli/Ozxw2PimUdHfOzMxEZWUlNm7ciLy8PLBYLPT19eG7777DvvXWpcR++tp5AMBtm3JRm0bvZNCaJ8/Suj9n4ogIkZGQr+Nza0KAYyIkk8nQ3NwMgiDsGqpnDa4QIYIgIBAIMDQ0hNLSUqN2G8aOi4lMPIuqR49DacOp0jxp2QW7fcyaoQ0Al8tFbGwsVeIrk8kwPz+PX5YM46U289EQf/z7qqtXdyyPEfifd5txuM+8o4Ov4UrLHm+FESEdZmZm0NLSgoSEBBQVFTlcGmkKZ7s6KJVK8Pl8yOVy1NXVWR3SM5GQ50AQhE1Nn7YgUdiXCg4KCkJQUBB+sy0ZL7UZFioQ0C8NBy59rgFf/7ae+u9nf1yh9/eusTlc/3o77OFPn7XjviuL7XquK3EkEloNU1UBN6fjnNkrZMtFnhy93dzcjKKiIpSUlDhNgADnRkILCwtoaGgAl8u1SYDI42IiIc/AWQIEACNWjvQ2xze/rTZ4xPC3TGBySY2BgQEsLi4aPa8KU2LsrqT7on3Orue5GkcKE6RSKRMJeSu2iJBCoQCfz4dCobD5om0vzhKhsbExdHR02GwjRMJEQp7BJ+eNr9l4EglhwfDnAEqTP7Plc++at0ewp3QEA2LgorxgpKckIzo6Wq+wwJbheCRbi2xvEncHjq4JrYZIaFWLkFAoBJ/PR3R0NCorK11mW0O3CGm1WnR1dWFiYgJr166lcvi2wlSreQZ7vhxx6PlxwcCMlKaDMcO53dYZlj78/we9vj0gBdALY6k7W7n/6hKHnu8qHF0Tsve37E2synQcQRDo6+vD+fPnkZOTg7KyMpf6ptEpQnK5HGfPnoVIJEJdXZ1DJ60r03FarRaTk5MQiURM9KXD1mctz/GxxAOXZljcho41STabjTB/e37DK1N3thBvrTGdB+BoOo6JhLwUcyKkVCrR2toKsVhssmnT2dDVrCoUCtHc3IzY2FgUFxc7vI7lqnScSqVCc3MzxGIxtFotCIKgelFiYmLstkPydnZ/0IExsePnxW8+HrK4zT/PjGBHfabDr3XyDxfQsH5l6HVnWtgSQv3wzf/YNrvInWg0Grvnpq2WwoRVJUIikQjNzc0IDw9HfX2922xAHLXtIQgCQ0NDEAgEKCgoQFpaGi1RpSsiIbFYjKamJoSEhGDdunVgsViQSCSYm5vD+Pg4uru7ERISQglSRESE2/3DXBEdvnpyEF90zjr9dUj+fGSYFhGiP5thLEpafiw9wh+f3GF/Z747cLQ6jilMcDLOTMcpld/3MRAEgeHhYfT09CA3NxeZmZludUV2JB1Hjl8QCoWorq5GVFQUbcfl7EhoenoaLS0tSE9PR25uLtRqNbRaLcLDwxEeHo6srCyoVCpqxDU5vC0qKgoxMTF6NjK+hGBagmeODlu1bXwgMG25FcilnPnDOqx74oyT9r78O00OYWFPpRbHjx/XG+YXGhrq0Q7njqwJicViqwfaeTM+Hwmp1Wq0tbVhfn6e9ou2vdgrQlKpFDweD1wuF/X19bSnrZxVmEAQBPr7+9Hf34+SkhJqMJsx/Pz8kJCQgISEBGp429zcHGUjExQURAlSZGSkU0vpXcVNr/Gs2i4hEPjsznWopumCf25gGjVZjo9wDw4IQCgXEDtpjNeaxGD86xfVIAgCMpmM8robHBwEm83W87rztFSuo2tCzmiW9zR8WoSWlpbA4/EQFBSEDRs22J2bpRt7RIhspE1OTqbVx87wuOgWIfImQCQSYd26ddQ8emsgh7eFhYUhMzMTarUa8/PzmJubQ1dXF1QqFaKioqjUnTf+YPd/2g2Z2vK5kBzuh0O/WV4L4WB5/o6j/PztLrTucVyEAODUHzeh/OHjoDuO3pwZhL/etNyTxGKxEBwcjODgYKSmpkKr1WJxcZEyX+3s7KRSuZ5yk2JvOo4gCEgkEiYScjbOTMdJJBKcPn0aWVlZyMnJ8aiQnRQhgiAsHhdZyTcwMIDi4mIkJyc77bjoTsfJZDI0NTVRjbOO3qVyuVzExcUhLi4OBEFAKpVibm4Os7Oz6O3tpcw2Y2JiEBUV5fYLkCWkShX+zZ+yuF1ebCA+uO17F+RmO/pqXAGf5uO6o5SDa2tSTf6dzWYjMjISkZGRyM7OhkqlWnGTEhER4dbUHWPbYxmfi4Q0Gg3Gx8chlUpRVVXlkXX25ElpSYRUKhVaWlogFottjiLsgc4RE2TlHmmBRHfkxmKxEBISgpCQEKSnp0OtVkMkEmFubg49PT1QKBSIjIykxhZYOxzMlfz8Hy0WtwkL4OgJkKfTumcTHv28A2/z9Issnv+vHGwsTAEAyFRq1B5oMLsf3h/rwefzbfrO/Pz8EB8fj/j4eOomRTd1x+FwqMjZVak7R9JxTCTkhYjFYsp8NDg42CMFCAB1Upo7QclUYnBwsMsq+eiKhIaHh9Hd3Y2CggK9QVrORNdsk1w7mJubg1AoRH9/P/z8/ChBioqKcmlfmDFUKhXaJy375/zHRDXYs9vy8D8fC6x4JfMX8QufOI6jf3B8CJ0u91yxBvdcYfrvQX5ctNy70Whp91Pb83Fp8bLZriOTVXVvUtLS0qDVarGwsAChUIjR0VF0dnYiNDSUEqSIiAinRM72puOUSiXUajVTou1s6LwznZiYQFtbG9LT0xEVFYWuri7a9k035A/L1AWffC+ZmZnIzc112R28o4UJ5AiMyclJVFVVITraPdYqumsHaWlp0Gg0EIlEEAqF6Ovrg0wmo9I0MTExbknT/PLtNovb3F6fgpgQ4+uYF5ckAVaJkHnmFA7vwi5YLJZF3zg6x3uz2WxERUUhKioKOTk5UCqVmJ+fh1AoRGdnJ1QqlV7VnbXD5CxhrwiJxWIAYETIGyAta8bHx1FeXo74+HjMz8+7fHKpLZAnt6EIkUaqo6Oj1HtxJY4UJiiVSvB4PKjVaqsm0Lryos/hcKgoKC8vj6qwmpubw9DQEDgcDiVI0dHRRqNOuo+3cdT8OAUOgF9flGN2m8evzsQf/zPo8LGo1Gr4uTkyNAZBEE7rEfP399erwtRN3fX394PL5epV3dlT1EQQhN1rQmKxmLqZ8nU878yzAalUCj6fv2L2D5fL9XgRMqyQI41UlUol6urq3LIgaW86bnFxEU1NTYiMjERVVZXVqS53rdEEBQUhJSUFKSkpVJpmbm4Og4OD6OjoQFhYGCVazsjJ84ctO0DfuzXb4jZXlKXTIkJ3vNuCV35a6fB+6IbOSMgc5lJ3IyMj6Ojo0EvdRUZGWnVc5G/JHhEiHbTd3ajtCrxWhKanp9Ha2oqkpCQUFhbqfVl0LrA7C10RIp0cIiMjXWqkauyYbI2EJicn0draiuzsbGRnZ3vc4r8ldNM0ubm5UCgUmJubw9zcHEZGRqhpo1qtFkqlkpYy/98d7La4zQ1VpqvCdLHHgXrFPkbEDj3fWbhKhAwxlrojo6SOjg6o1Wq9AgdTRS/kNcie9yAWiz2ymMYZeN2akFarhUAgwPDwMNX4aAiHw6FCYU+9kyB7mUZGRtDV1eURTg62rAkRBIHe3l4MDg66JXXoLAICApCcnIzk5GSqD2VychIEQeDEiRMroiR7zq9ZmfmuzspU29YBbt+QgBdPWi71NkWgZ7TPrUCr1XrERdjf3x+JiYlITEyk+nfIdG5fXx/8/PwoQYqKiqJuVBwRodVSng14QCRky4VPLpeDz+dDpVKZnf1DRhKOlEc6GxaLhb6+PiwsLKCyshIxMTHuPiSr03FqtZoqHV+/fr3PlpGSfSh+fn6YmppCXV0dVXHX0tLiNOPVl35SatP2v76wwCERCg/2LJcBEk+8iWSxWAgNDUVoaCjS09Oh0Wio1N3Q0BDa29sRFhZGRUhsNtsuIRWLxbQVR3g6bhcha5mbmwOfz0dsbKzFdQfdEmh3mZSaQyaTQalUQiqVoq6uzmP80KxJx0kkEvB4PAQEBGD9+vUOpae87Qfm7++PpKQkynZoaWmJMl7t6uqi1g3MGa8OzpovSACAABefs6H+nnkZcGZhAl2QRS1kJahCoaCq7sbGxqDVaqmZZeZSd4aslqmqgBeIkK5jQFFREVJSUix+ieTdhyeuC83NzaG5uRkcDgd5eXkeI0CA5UhodnYWfD6fNusgbx6gR64V2Wq82jVJ//qLROmYaVusB87nIR1FPF2EDAkICKBSd+QaUlRUFGZnZ1ek7kxVYgLfR0KrAbd/w+YERalU4vz58xgfH8e6deuQmppq9d2zp1XIEQSBgYEBNDU1IT8/3yNLL02lRgmCwODgIHg8HgoLC53igODtkMara9aswYYNG1BZWYnw8HBMTU3h9OnTOH36NAQCAfIiLAuvVK60uI0ujQPz9h42ACAjyvPSceR56M3nmVarhZ+fH9LT07F27Vps3LgRhYWF4HK5GBwcxHfffYdz586hr69vxXBHe2cJHThwACwWCywWC6dPnza6zeLiInbt2oWMjAwEBAQgIyMDu3btwuLiosn9vv3226itrUVISAiioqJwxRVXoLGx0ebjM4bHRkLz8/Pg8/mIiIhAXV2dzWk1T6qQ0zXxJAfpjY+Pe9xEUWPpOI1Gg/b2dszNzaGmpgaRkZHuOTgvwpjxKlldtfug5fHWhI0X3i/a7V8PAoCKNNcPdrQE+dvwdhHSLc/W7VcDllN35HnR2tpKbX/y5ElMTU3ZfKPa2dmJBx54ACEhIZBIjLtxSCQSbN68Gc3Nzbj00ktx4403gs/n4+mnn8aRI0dw4sSJFRHYI488gj179iA9PR233347xGIx3n33XWzYsAGHDh3ChRdeaNsHY4DHiRB5193b24u8vDxkZGTYtXZgacS3qyDXUPz9/VFfX0+toXiSSJIYpuPkcjl4vOUxA3V1dQgMDHTXoXk1XC4X8fHxUHKD0SudhqVJoiE2rtGcH7G8zmSO6kz3OFuYgzwPvW3dUBdLhVEBAQF6a4xisRjnz5/Hl19+SVl23XbbbdiyZQsuvvhis2NoNBoNfvazn6G8vBz5+fl46623jG534MABNDc3Y/fu3Xj88cepx/fu3Yv9+/fjwIED2LdvH/W4QCDA3r17kZ+fj7Nnz1KTqH/729+itrYWO3fuRFdXl0NtJR51m0GOfR4aGkJNTY1DJcuekI6bnp7GqVOnEBcXh+rqar1FfEenqzoD3UhIJBLh1KlTCAkJQW1tLSNANHD5C01GHjU8v7Xg8XgYHh6GRCKxat1sWqxy6Lgigz2vRtsXIiFbLHvI6PnCCy/EkSNH8N///d+49NJLERISgr179yInJ8fs9ezxxx8Hn8/H3//+d5OvSRAEXnnlFYSGhuKBBx7Q+9s999yDqKgovPrqq3rn3GuvvQa1Wo09e/ZQAgQAxcXF2LFjB/r6+nD48GGr3qMp3P4NkyKzuLiIhoYGaDQa1NfXO5z2cWekQRAEBAIB+Hw+SkpKjC7iOzJd1VmQkdDY2BjOnTuHzMxMlJaWevxIBG/goybrJqfuuzQFsbGxmJ+fx7lz59DQ0ICuri7MzMxArTZegOBZZxE9eGJ5tq04MtpbqVSipKQETz31FNra2tDX12dyX21tbdi3bx/uu+8+FBcXm9ynQCDA+Pg4NmzYsCLlFhgYiE2bNmFsbAy9vb3U40ePHgUAbNmyZcX+LrvsMgDAsWPHbH17erg9HUcQBNWwSWfXvbvScUqlEi0tLZBKpWZ7aDxRhIDl4+/s7MTatWs91oXcG7n/i0GL2wRzgR/W5gGAWeNVsuLO00dbO4I3VsYZ4sgsIcOpqqZScWq1GjfffDOKiopw9913m92nQLBseJuXl2f07+TjAoFA79+hoaFITEw0u70juF2EyK57uhs23ZGOW1xcBI/HQ1hYmMViCk8TIZVKBYFAAI1GgwsuuMAl5aHecgF1tJT8y7ZJq7Y780d9V2lzxqvkfBxPaHJ2Br4SCTkyS8ia6rhHHnkEfD4fZ86csVi8tbCwAAB6aTVdyHll5Hbkv025oRjb3h7cLkKpqalITk6mfc3B1em48fFxtLe3Wx3NeVJhglgsRlNTEwICAsDlcldNf4Kr+N+Peyxu89qNhRa3MTReFYlEmJ0lh8eZr7gzhb+HZlo9xbLHERxJx1kjQnw+H3/605/whz/8AZWVnmdAay1uv9UIDg52yqK3q9JxWq0WHR0d6OzsREVFhdWjxD2lMGF6ehqnT59GYmIiCgstXwgZbGN2UWpxG382UJ1tm/cem81GdHQ0pIFx//8R+y7Y15R4puefr0RCzhShn/3sZ8jJycGDDz5o1T7JCMhU5EL2CelGShERETZtbw9uj4SchStEiPSyU6vVqKurs6mun81mQ6VyrKrJEQiCQH9/P/r7+ykj2KWlJZcKo7ff6VrDD/9muaHvs19X273/r7pmrNjKdJR03xX5dr+2M/EFEbJ3TYicb2QpI8Hn8wHA5E18XV0dAODDDz/E9u3bLa7hGFszysvLw6lTpzA5ObliXcjSGpO1uF2EnHUh4nK5UCicNzZyfn4ezc3NiImJQXFxsc0nmzvXhDQaDVpbWyESibBu3Toqt+voZFWGlcxbYX6QGO6Ie4Y1vx/T26hVKnBoMl6lE18oTHBkTUgsFls0Bv7FL35h9PHjx49DIBDgmmuuQVxcHDIzMwEsi0VycjJOnjy5wqVbLpfj+PHjSE5ORm5uLvX45s2bcerUKXz11VfYsWOH3uscOnSI2sYR3C5CzsJZay4EQWB4eBg9PT3Iz89Henq6XULqLhGSyWTg8XjgcDioq6vTc352ZLIqw0oaeqYtbvOX6woceo1ryhLxj3MTdj//w0MnkZ1o2XjV1fhCJOTsdNwrr7xi9PGbb74ZAoEA99xzD9avX089zmKxsHPnTuzfvx/79+/Xa1Z99NFHMT8/j9/85jd617NbbrkFTzzxBB5++GFs27aNSr21t7fjzTffRE5ODi6++GK73iOJz4qQM9JxuhY21dXVZjuYLeEOERIKhWhubkZCQoJR/zcyEiIIwiWpMvK1fJXb/t1lcZsLCxIceo3CRMfGaFRVlyOcozZqvBoTE+O2JuXVXJhApuPs8Y6zxO7du/HJJ5/gwIED4PF4qKqqAp/PxxdffIGKigrs3r1bb/v8/Hw8+OCDuO+++1BWVobrrrsOEokE77zzDlQqFV5++WWHh3C6XYScmY6jU4SkUimam5vBZrNpsbBxdXXc8PAwuru7UVBQgPT0dJPHBMBlIuRNOOPzCPOALFhe8nKJd0JCAmUdMzc3h8nJSfT09CA4OFgvSnJV47IvREL2rgnJ5XJoNBqniFBISAiOHj2Kffv24eDBgzh69CgSExNx1113Ye/evUbXofbs2YPMzEw888wzeOGFFygLsv3796OmpsbhY3K7CDkLOi/y5AgDY6PE7cVV1XFarRadnZ2YnJxEVVUVNffEGOSF1hcuAN7AFYWRtOznnztKcNObbQ7vx5Tx6tzcHDo7O6FSqfTGUzjTCd4XzkF714RI81F7Rej111/H66+/bvLvEREReOqpp/DUU09Zvc+bbroJN910k13HYwmfFSE6IiHdCrI1a9YgJSWFpqNzTTpOqVSCx+NZXb2nGwk5G4IgoNFooNFoqMiLxWJ5/YXHFv7FF+G+qxzfT1laNFIjAjC6YFshzpm7zFflkcar8fHxemOtZ2ZmIBAIEBgYSAlSVFQUrVGSrxQm2POZiMVisNlsj5o15kzcLkLOSvs4uiZEjrBeWlrSqyCjC2eL0OLiIpqamhAZGWlxEi0J+V04W4QIgtAzqNT9b1KIVpsgOcoXd67DJc80YEpi3YC75rs32HSBNBxrrVarIRKJMDc3h56eHiiVSkRGRlKpO2sniJrCFyIhe9NxZHn2akmJu12EnIUj6TixWAwej4egoCDU1dU5NMLaFM4UocnJSbS2ttrsxaebjnMWZARE/kC5XC41RZN8XPd7Y7PZ1P98Dbp/fN/8rh7HBdP49XumCyLuuCANv9qc5fBrcblcxMbGIjY2FgRBQCaTYW5uDnNzc+jv74e/vz8lSFFRUTYvXnt7YQJ5Pttz3pJTVb35/duCR4iQM/pT7E3HkRfwjIwM5OXlOe1EcEZhAkEQlBdfWVkZEhJsq7xydiSk1Wqp/5HRDvB9GpC8aySFiBQn0j3am6Ika77bb++soP11N+XFo3XPsgvCZw2t4E1KsGlNDjbkRTutqIDFYiE4OBjBwcF6xqtzc3MrjFdjYmKsusB6eyRE3sjZ85kb9vD4Oh4hQs7A1nScVquFQCDAyMiIXRdwW+FwOLRe7Mn0oVgsNuvebQ5yXYbuSIgsxSa/D10BMoZu5KMrXN4UJX3SPGpxm+gIelO8huTFBqAgPgi5uXGWN6YRwwmiZJQkFAr1jFfJKMmY8aa3ixB5njoiQkwk5OXorjVYOpmVSiX4fD7kcjnWr1/vlNJIY8dHVyQklUopA9L169c7lD6ku2GVFB9yn6TQ2XI8uoJEfqeeHiU98OWQuw/BY0rtg4KCkJqaitTUVMp4VSgUYmBgAO3t7QgPD6dSd2FhYVRmxJvnWOnecNkKEwm5AWel4wDLZZILCwvg8XiIiIhAXV2dw41X1kLXmhBZPp6cnGx0eJ6t0Pld6K7zWIp+rMFU2o4UJm+JklyFp4iQLqTxanR0NHJzcyGXy6ky8OHhYervCoWC9mIgV0Kuedrz+YvFYpfcCHsKHiFCzoC8+Gg0GpNzNkZHR9HZ2Ync3FyHRonbe3zk3bw9r0sQBIaGhiAQCFBUVITU1FRajouudBzdAmQMY1GSrii5q+Lu5eMDFrfJinZ+p6onipAhgYGBSE5ORnJyMrRaLRYXF6kCB5FIhPn5eSp1Fx4e7vHvh8TRWUJMJOQDkBcdYykv3QZOuofp2XJ8gH0XCq1Wi/b2dszOzqKmpsbhUeiGx+VoJGSqAMGZWBsl6fYjOUuQnvtuxOI2b/13qVNeWxdvECFd2Gw2IiMjERkZSdnWBAYGQigUUo7RZNouOjpaz/fQ03DEN84aB21fwiNEyJXWPXK5HDweDwRBoL6+3m0NYbp38LZcDMnjB0CLfZAhjkRCthYgOBNTUZKuQOoeo6vXksJDnec2QOLNvnxarRZ+fn5ISkpCUlISCIKgoqSxsTF0dnYiNDSUEiRPMV4lcWS0N5OO8yEMK+RIA8+4uDisWbPGrQuf5GtrNBqr16FEIhF4PJ7d4yOswd5IyFgKzFPuwi2VgNOZtnvtO8upuA9vXmPXvm3Fm10HDG/OWCwWIiIiEBERgezsbKhUKqrirq2tDVqtllprcqfxKomj6ThHzJG9DZ8WIbICTXf9pLCwEKmpqW6/QNraGDo2NoaOjg6nr1/ZEwnpRhmeUp1mDnNRkrHiBlt45rjlVFxuSqxtB2wn3h4Jmfvs/fz8kJiYiMTERIvGq5GRkS4/Jx1Jx8lkMqSlpdF8RJ6LT4sQl8ulyq/n5+dpXz9xBPJibemCr9Vq0dPTg9HRUVRUVCAuzrk9H7ZWx7miAMGZ6EZJ5Hdh2CirVCqpxy2JrCX5rk1zbPSCLXjbmpAutkRxhsarKpUK8/PzbjFeJXFEhMRisUuO0VPwCBFy5g+lp6cHQUFBqK+v97iFTEsipFKp0NzcDLlcjrq6OpcsVtqSjvN2ATKEvOjpRkkymQy9vb0IDw+3WAL+bYfl4XIv3lhC81GbxptFyJFmVT8/vxXGq3Nzcy4xXiVxZE3ImoF2voRHiJAzmJ6extLSEqKjo1FVVeWRKSJzIiQWi9HU1ISQkBCX9i9Zk44jCxB0K+AWFSr8u3EMLx4ZhFRHw1gAPrx1LfKTIqnHRFIlWscWoNVokJEQhoxIx8wunYVUKgWPx0NUVBQKCwupz8ZUo+zvPhRY3KepdgFn4O0iRMex6xqvZmRkQK1WY35+HkKh0CnGqySOrgkxIuTFEASBvr4+DAwMIDw8HHFxcR4pQIBpEZqenkZLSwvS09Od6l9nDEvpOMMChPv/04WPWkyPsSYAbH+Zh4gANhYUpsWNywIOXFuErSWJdh87nZBlwWlpacjJyaG+A2saZU2RGUW/Ea45vF2EnPG75XK5iIuLQ1xcHDXBlGyWpcN4lcTREm1GhFwMXT8UlUqFlpYWSCQSrF+/HgMDAy6dXmorhtV7uvOLSkpKkJSU5PJjMhed6TaDvn5qCE98O2j1fs0JEACoCWDXB534x5kR/PMXjk9rdISJiQl0dHSgoKDAZBOwYXHD++ctFyT845a1tB6nJRgRMg+LxUJISAhCQkJoM17VPX57R3szkZCXsrS0BB6PR6Wv/Pz8HJ4p5Gx0L/gajQatra0QiUSora1FRESEW47JVCSku/5z5f+dwbDItgFq1sIbE+P1U0O4uS7DKfs3B0EQGBwcxMDAAMrLyxEba10VG5vNxt4v+i1uF8JlW1XcQCfeKkLuKC+nw3iVxJxTiyUYEfJCxsfH0d7ejqysLL3UiTPGJdAJKUIymQw8Hg8cDgd1dXVuLaAwVphACpBSqcS6PzdA6eSp5H/+ut/lIqTVatHV1YWZmRnU1NTY5UJuDvKe2FQJuDMuuEwk5Bj2GK+SMGtC1uMRImTvD0Wr1aK7uxtjY2MoLy9HfHy83t+5XC4UCufcsdMBm83G0tISOjo6EB8fjzVr1rj9h6dbmKBbgCBVKLD+wCm4QtJd3d2iVqvR2toKmUyG2tpap7hofPjLSgQGBrp0VpK3i5AnHbu1xqtk1Z29a0JkOo6x7fECFAoFmpuboVKpUF9fb7Su3tPTcUqlEv39/SgqKkJ6erq7DwfA9+k43QIElUqFOhcJEIlSo4G/CxwtFAoFeDwe/Pz8UFNTY1cKZWJ21uI2uYnL6VVXzkryVhEizz9335CZw5Tx6vDwMDo6OigH7bCwMJuMV6VSKQiCoD0S92S8UoRI+xqy/NpUBYunpuPI1I9UKkVaWppdAvRt1wT2/UeAObESWgBVaaF47WfV8PNz7CslU4S6M4Au/es5qB3aq+182TKCa9ZmOvU1yDHuUVFRdkeh//3ySZyfUNr1+s6eleStIqTr6+cN6Bqv5uTkQKFQ4Pz581CpVDYbr0okEgBg0nGuxtofCkEQGBkZQXd3N/Ly8pCRkWH2ufaO+HYmSqUSPB4ParUaMTExNqd+TvfP4edv8qAxyFmdGxGj5E9HsT49HG/8otbu49PthWGz2fjZ6+cxK1HZvT97ufs/A/iAP4vXb652yv7JEuz09HRkZ2fbdbFes/8IbcdjzN/O0SjJW0WIvPnxFhEyJCAgABwOBxkZGYiJibHJeFUikYDL5XpcY70z8ZpvWaPRoK2tDb29vaiqqrLKP83T0nGLi4toaGiAv78/1q1bB39/f5t82o50TeNnb6wUIF1ODy/i0ifsuzgSBAE/Pz8MDg6ira0Nz3/VgsaRJaueW5pgXkxLEwLx7s2l6HjgInQ8cBFKkiznvM8OL2HN/iNY99gRHOs23YtkKxMTE+DxeMjPz9crZLEFOgXIGGw2G1wuF/7+/tT/OBwOVThC2gmpVCpqjckQb/WO87ZIyBjkmhBpvJqdnY3q6mpccMEFSE9Ph0KhQFtbG7777ju0trZibGwM09PTkEgkCA4OtvjeRSIRfvvb36Kurg6JiYkICAhASkoKLr74Yrz//vtGv/vFxUXs2rULGRkZCAgIQEZGBnbt2oXFxUWTr/P222+jtrYWISEhiIqKwhVXXIHGxkaHPx9dPCISAsw3SZLVYywWC/X19VY75HqSCE1OTqK1tRXZ2dnUnbct01XX7fsGIiv1anhJg3+eHsRN6zOtPj7yjjs7OxsJCQnoGJzE/50Zx3KZgPmL9P9cmI7bNuVAo9HgF/9oxtnh5ZM6OpiDN35aipzElY7AB64twhXPW3cyLymBO/7VDqAdgRzgrosz8dP1tpu4kiXYg4ODNpVg66LRaFD68HGbnvPuzY7NDrJ3oqy3RkK6rubeiqnCBH9/f6PGqxMTE/iv//ovBAYGQq1W49tvv8XGjRtNRkSzs7P4+9//jvXr12P79u2Ijo7G9PQ0/vOf/+C6667DrbfeipdeeonaXiKRYPPmzWhubsall16KG2+8EXw+H08//TSOHDmCEydOrCiGeOSRR7Bnzx6kp6fj9ttvh1gsxrvvvosNGzbg0KFDuPDCC2n5rFiEh9wuKZVKoyJEjq9OTExEUVGRTXdHc3NzaG9vx6ZNm+g8VJsgCAK9vb0YHBxEWVkZEhISqL+RC5gFBQVm91Gw9xu7Xrt73yVWHZ/uDCByBEPtY8cgXlGLvVKQdl2ciZ0XZNl1fBV/OuJQubcfG9iYE4m7fpCNnHjzfVXkOtzs7CzWrl1r18Ivb3gWN73eavPzOh64yObnWIuhC7jub6ixsREFBQWIioryqqhCKpXizJkzuOgi531uzubYsWOorq62qcptdnYWTzzxBN566y2EhYVBJBLh4osvxtNPP42cnBy9bcmbEMP18KWlJaxfvx4dHR1oa2tDcXExAGDv3r3Yv38/du/ejccff5zannz8gQcewL59+6jHBQIB1qxZg+zsbJw9e5bqW2xvb0dtbS2SkpLQ1dVFi52Yx56ZpHsAj8dDQUEBiouLbf4huTsSUqvV4PF4mJiYwPr16/UECLBsYArYL0AAcOBQu9m/kxcvwyF0f/q824gAAYYCtD5Oi+rQBYyOjtpVCv/lb9bZ/BxdVFrgsECEq19swpr9R1D/5+/wxFfdWJDK9bZTq9Vobm7GwsICamtr7RKgdY8c8TgBApa/Mw6HA39/fwQGBiIgIABcLhfT09NQKpXgcrlQq9Vm03aehif0CDkC+buy9T3Exsaivr4eWVlZGBkZwalTp7Bx40ZER0ev2JbD4RgVgLCwMFx22WUAgN7eXup4XnnlFYSGhuKBBx7Q2/6ee+5BVFQUXn31Vb0bmNdeew1qtRp79uzRa5wvLi7Gjh070NfXh8OHD9v0/kzhMd+0buhNXjSGh4dRW1tr0jrFEu4UIalUitOnT0Oj0WD9+vVGL3zmjk+pVDokQADwaoNpV2fdSizdqiuVSoV3Gsct7rsgLgjP/nTZ2WF8fBzfffcdzpw5g76+PiwuLlq1HpEYEYynrs23/g1ZQCRT4++nx1H3xCms2X8EG5/8Dk8e6sap06dBEASqq6ttHnb2ZesE1uw/giUbywMT/Z0vQMZgs9kYGxtDf38/ysvLERERoTdAkVxLUqvVKyInT8HTy7MtQQq9PX1C5FRVFouF0tJS/O///q9NA+7kcjkOHz4MFouFNWuWhycKBAKMj49jw4YNKyKzwMBAbNq0CWNjY5RoAcDRo0cBAFu2bFnxGqTIHTt2zNa3ZxSPWRMiIctmAwMDUV9fD39/+00f3SVCZAoxOTkZBQUFJn9Q5iIhW9cdTNExIsSaNP07KXMjGG59u9Vis2ggl4X3b6sFm81GWFgYsrKyoFQqMTs7i5mZGQwNDYHL5SI2NhZxcXGIjo42+YPcWpqC7NgQbH+Z5+hbXcGcRI1Xz4zjVQCACpGHT+Pq0lj8ckMWYsLMi9GcWIoLnzpjV2/Ut78qQ1JsjB3PdAyCICAQCDAxMYGqqiqEh4cDMD1R1pmNso7g7ZGQIyJka6OqSCTCM888A61Wi+npaXz++ecYGRnB3r17kZeXB2BZhABQ/22I7na6/w4NDUVi4kpDYcP9OopHidDU1BRaW1uRlpaG/Px8hxcmORwOdcfvipNad4JrUVGRxQjOlAipVPSVRP/s9Sacu//7tSHdsl9DAdJqtTg7tGBxn/+5o2bF5+nv76/XvDc/P4/Z2Vl0d3dDoVAgOjqaEiXDaCQ/KRLt91+IO97h43jvvIPv2DQimRr/ODuJf5ydBAAEcIDqzEg8dFUhEiOWq/tUKhXKHz1h92ucuavaLY2GWq0WnZ2dEAqFJtciXNko6wjeLkK66W1bsdWyRyQS6a3l+Pn54c9//jN+//vfU48tLCz/pk35UZI3K+R25L8NHWjMbe8IHiNCvb296OvrQ2lpqVH1tQfdNISzT2qtVov29nbMzs5aPcHVlAg9/WUbbce1+P/TSIYFCMaG0B3unrG4v19tTEdKlPk7NTabTRk95ufnQyKRYHZ2FpOTk+ju7kZoaCglSGQ3OYvFwos/qQAA3P9hC95vnbP9zephuapPoQFO9olw8bOnHXytZdyRfgOWz++WlhbI5XLU1NRYlXJ0dqOsI/iCCNk75NHWSCgzM5P6XY+MjODdd9/Fnj170NDQgPfee89lc8gcwWOOMDw8HHV1dbR2CuuKkDOHicnlcvB4y+mkuro6q9cdjImQQqXCq43WXYBfvT4Vez8fw6jEfAJNdwQD8H0FnCHPHh40u5+oIC7uvCjH7DaG6A4VI0cvk2m7pqYmsNlsvbQdl8vFQ9eW4aFrgUWpDLs/6saJ3nmLY7ONvLLNz3CEz+6odOnrkZDTdwGgurrarvPc3hJwZ+FpvnG24sgsIXvNSzkcDjIzM3H33XeDw+Fg9+7dePnll3HHHXdQEZCpyIXsE9KNlCIiImza3hE8RoQSEhJoX78h796cuS60sLCApqYmxMTEoLi42KaTz9ixbXvhrFXPbfzDOoSFhWHs36MWtzWVfjNkWCQ1u58Xf1Ju1bGZw8/PD0lJSUhKSqKciWdnZyEQCCCTyfTSduHBQVR0RBAE3mscwd8ahjC54GoTIfNszI5CVpzrR28oFAo0NTUhMDAQZWVltI2pNhYl6YqSbh+PM6Ikby9McESExGKxXf1rumzZsgW7d+/G0aNHcccdd1hcwzG2ZpSXl4dTp05hcnJyRWbK0hqTrXiMCDnrzseZxQnkCInc3FyrHByMHZtuJPQhbxQDczKLz/vLtZnUuoM1tU2tw7MoSYuxeHwqCx9TaUq4Fa9mPbrOxLppu+npafT09CAkJIQSpIiICPyoJh0/qln22ZOpVPjnmTH84+woZsSutxUiee66IlyyxvXTYKVSKZqamhAZGelU93VroyQyuqYjSnJF+tyZ2DvQDqBnqur4+HJ1K5mKy8vLQ3JyMk6ePLki3SeXy3H8+HEkJycjNzeXenzz5s04deoUvvrqK+zYsUNv/4cOHaK2oQPv/aatxBkiRBAEurq60NHRgYqKCmRlZdklorrpOLVajbs/6rL4nEg/YEtFrsXtdPnR661Q0ljs4CxCQkKQkZGB6upqbN68GVlZWZRb+rFjx9DW1oapqSmoVCoE+flh5wWZOPSravx9SxBeuCQQd1yQirhg19xXRbCW14DcIUCLi4s4d+4c4uPj7eqfcwQ2mw0/P78VdkLA8sVXtwTc3r4kX4iE7D1+qVRq1ZoQ2fdmiFAoxL333gsAuPzyywEs3yDs3LkTYrEY+/fv19v+0Ucfxfz8PHbu3Kl3DbvlllvA5XLx8MMP671Oe3s73nzzTeTk5ODiiy+26z0a4jGRkLOgW4RIZ1yZTIa6ujqH5n7oitClz1hXkfXdbvvuPrb85QyO/X6jXc91B35+fnr2JmTarq+vD62trYiKikJ4eDjGx8cRGxuL2qIibGaz8ZuLl1MEI0Ip3jk3gk9apyCU0vP931yTiBvWpSEz2n0Ox6TxalZWFjIzM912HIDpKEm30AGwPW3nC4UJjqTjrImEXn/9dbzyyiu46KKLkJGRgZCQEAwNDeGzzz6DWCzGf/3Xf+EnP/kJtf3u3bvxySef4MCBA+DxeKiqqgKfz8cXX3yBiooK7N69W2//+fn5ePDBB3HfffehrKwM1113HSQSCd555x2oVCq8/PLLtBU9eIwIeUM6TiwWo6mpCSEhIVi/fr3FReAlmRwBXA78TWxHitBxwRTGreiGfO66Ivj76++rIjkUzeNii8+dkai91kuMxWIhKioKUVFRyMvLg1QqxdDQEIaGhiiB6u3tpdJ2bDYbadHB2H1ZAXZftmyJNDizhL+fGcHnrTOQqmy7O48JAL77o/stZKanp9HW1oaCggKkpKS4+3BWYGotydbiBm8vTHAkHWdtYcJ1112HhYUFnD59GsePH4dUKkV0dDQuuOAC7NixAz/+8Y/1PsOQkBAcPXoU+/btw8GDB3H06FEkJibirrvuwt69e43eTO/ZsweZmZl45pln8MILL8Df3x/19fXYv38/ampq7Hp/xvAY7zhyeBrdnDlzBmlpaUhOTnZoP9PT02hpaUF6ejry8vLM/kiMOR0UxXHw0Z36FzJyLtIdRzUW13bWxAfjw1/Xr3icIAgUPvitVe9h35V5uL7KdO9S+Z+OwNz1+W83lmBjXpxVr+VMxsfH0dnZiaKiIsTHx2Nubo6quCMIArGxsdT/TN0oTC/J8KMXT2PKzBJcmD9w5m73iw8AjI6OoqenByUlJSb7NzwVMiLSjZJ0KzUNo6TBwUFIJBLK98zbGB4exsLCAkpLbTeuraurw0MPPYRrr73WCUfmmXhMJOQsHI2ESA+7/v5+lJSUICkpyeS2izIFah77zujfOmc02PzoNzh2z/eNoxwOB+90WBYgAEYFCFj+ER/YlovdH/ca/bsuL58YMCtCcaH+GF80PaDtwNf9bhUhgiAwMDCAoaEhVFRUICZm2ZUgISEBCQkJIAgCi4uLmJmZweDgINrb2xEZGUkVNwQHB1M3D/FhQTjyv54hMObQdf5eu3atTRYungIpLtY2ynp7YYKjo71X01RVgBEhs2g0GrS2tkIkEqG2ttZsXXzP+Dyu/tt5s/ub1PfVBIvFwrFpLSzVh7x181qTfyMIAleWpWJzThTWPXXO/OtbKG2+siQRLzcMm/x736z5Em5nQjoCzM3NobrauCsBObslIiICubm5kMlkVITU19eHgIAAxMXFIS4uDpGRkR5/oSMIAj09PZicnDT5nr0Rc42yWq0WUqkUXC4XKpXKo+yErMUdfULejMd8s562JiSTyXDmzBkoFArU1dWZFaBPeCMWBYjktn98LxRfdE7DUlNlVnQQarJW+pCRP1y1enmtJzQkBP/eWWF2X5ZWnX57UaaFLYDGQefZ6piCNLRdXFy0yQU7KCgIaWlpqKysxIUXXoj8/HzqxuLYsWNoaWnB+Pg4lEr7xnM7E9KBY2ZmBjU1NT4jQIaQLuBkxd3IyAhEIhG15mXMdNXTcXeJtrfBREJGmJ+fB4/HQ3x8vMUejIc+68RbZ8es3veJ/u/LHR/6rN/i9v/51cpR3caaBlksFvLiHDt5ORwO/AGYuyT/5l8tOPVHevoDrIF0owgICEBNTY3dFTkcDgfx8fGIj4+n0nazs7MYHh5GR0cHIiIiqLRdSEiIWxfGDW14VsOoZ4Ig0N3djenpadTU1CAkJGRFcYPuOU9GSJ4YJWk0GrvOU61Wa7Ntjy/AiJABIyMj6OrqQkFBAdLS0sxejG546Qz4Y9aNvyYJDfj+DklhYTGoOiNsxcK6oQWP7g+QS0PH/PXVyfinmVEOCwot5EoVAv2dZ4NEsrS0BB6Ph5iYGJsHGppDN22Xk5MDuVxOpe36+/vh7+9Ppe1cPRBOpVJRU4TtteHxNgiCQEdHB+bn51FTU4OgoGUzWWeVgDsbe9NxEokEAHw26jWFx4iQM9Nx1gxcIydvTkxMoLKyklr0BoAvW8bw4vEBzEg1yIkJwINXFWDHGy2YkdhezffnH1pf8fPmz6r1/tvcCAa6+N8tOWZFCAB+/W4bXt1hep2KDubm5tDS0oKMjAy7m4GtJTAwEKmpqUhNTYVGo4FQKMTs7Cza29uhVqsRExNDVds5MyqRy+VoampCcHAwSktLabPh8WTItOPi4qLFeU+2loCT/3Y19qbjSBFiIiEfw5pISKlUgsfjQa1Wo66uDsHBwQCAP33WgX+c1b8gz0pUuPx569Z/DGED2JRvfXmt7onsCgECAH8uF+EBHCwqTH9mpwZFTnltEt0SbEdL622Fw+FQUVBhYSHEYjFmZmYwNjaGzs5OhIeHU2k7cvgYHUgkEjQ1NSE6OprWqM+T0Wq1aG1thVQqRXV1tU0CbylKcqcLuL2RkFQqhZ+f36pIv+riUSLEYrFon/RoSYQWFxfR1NSEiIgIVFVVUbncrc+dtMrHzRa++V2dzc8heyqsMiGl6YL47A0luOUffLPbjM1LkRIVTMvrkZDl8MPDw3ol2O6CxWIhLCwMYWFhyM7OhkKhwOzsLGZnZzE4OAgul4u4uDjExsaaHdxniYWFBfB4PKSkpCA3N9erGzWtRaPRgM/nQ6VSoaqqyqHhlYBnzUqyt8RcLBa7fT3SHXiUCDkDcyI0OTmJ1tZWZGVlIScnh/ry9/+nk3YBunJNrIk5PKZPOMMCBEsREJumk3ddVjRYMG+O+rt/8fHv220XVVNYU4LtbgICApCSkoKUlBRqcN/MzAy6urqgVCoRHR1NiZK14zzm5ubA5/ORk5ODjIwMJ78Dz4CsdtRqtaisrKR93cvds5LsjYRIEVptrEoRIggCvb29GBwcRFlZGRISEvT+/lXnNO3H8dSPKmx+jqkCBFdwRVEMPus0PdeoY1pu8m+2olarwefzoVQqUVtba/UF3J3oDu4rKCiARCLBzMwMJiYm0NXVhdDQUEqQyMF9hkxNTaGtrc0taUd3QRZecDgcVFZWOn3omrG0nbOjJEfScXSmeL0FjxIhV6Tj1Go1WlpasLS0hPXr1xu945ZYmmlgIznR9uV4nbH+82HTCK6tTLO43cPbCvFZ50mTfydAj9EkXSXY7kR3cF9WVhaUSiWVthseHgabzaYEKSYmBhwOByMjIxAIBCgrK0NcnPutkFyBUqlEU1MT/P39UV5e7pbCC3OzkuiKkhwpTCDXo1cT3veLtxFdESJnsAQEBKCurs5kHjomxB9SJX13+q/+zM5KMhMTUM3hxzE/F+i+T3utEiFrcvSLMiUiQ+yPWsgS7NjYWBQWFvrMYry/vz+Sk5ORnJxMDe6bmZlBT08PFAoFAgMDIZfLUVJSsmoEiBzAR1b+ecJ37ayJso6sCa22RlXAgxwTnAUpQnNzczh16hRiYmIsLoT+vD4dHBoj4qRI4yeWpdeYWrRdCK8tNz/fhgAgockh4KVjA3Y/d25uDo2NjUhNTfXpajBycF9BQQHq6+sRHx8PhUKB0NBQtLa24tSpU+jt7YVIJKI9C+ApyOVyNDY2IjQ01GMEyBimZiWxWCyrZyWRkZW9kdBqXBPyqLPBGblQNpsNlUqFpqYmFBQUWHXB+0ltOm7dkIEQf8c/Hj8zuwiwcJ42DQltfr37tloeufuz15qt2leQuYMH8HrjpFX7MWRsbAzNzc0oLCxEdnb2qsiB6/bD1NXVYd26ddi8eTMyMzMhk8nA4/Fw/PhxtLe3Y2pqikoLeTsymQyNjY2IjIxESUmJxwqQIaSdkL+/PwIDAxEQEAAulws2m00JjVqthkql0hMk8kaCESHr8el0nFarRV9fH7RaLdatW2eTA/Fdl+bhzgszUfKnYw4dw8560xVPQf5+kKpNN7z2zdheocflchEeyMGi3HROrmNKApVaDT8L6y958UFoGZOY3aZjfB5rkq37XHVLsNeuXYvo6GirnuftkOXISqUSNTU1VBTu5+eHpKQkJCUlQavVYmFhgTJbJQf3kT1LpIuANyGRSHD+/HnEx8ejoKDAq282rG2U1R13biur0bwU8LBIiE7kcjnOnDkDqXTZ+dmcAakpPuJPOXwcv7vUdGQSFWy+NHVYaN+61PM/tjzH5EevWG64fe76MovbXPdKszWHREUC4+PjqKmpWTUCpFQqcf78eRAEgerqapNpYDabjaioKOTn56O+vh719fWIjY3FzMwMTp48iYaGBggEAszPz3tF2k4sFqOxsRFJSUleL0CGmIuS5ufnqWGVhlGSJZhIyAOg60RdWFhAU1MTVT575MgRuxYLxxccK04wl80iCALBXMOLCQHdviGFnSMoKtOjLPb5dE1LoVSpTE59BYD4cOuKDna+2YhXdlSb/LtKpUJLSwtUKhVqamq8ogSbDkgbnpCQEJvXQoKDg5GRkYGMjAyoVCoIhULMzMyAz19uIiZthGJiYjzOX45sAE9LS1sV6VYySiKLT9asWQMul2tzcYNEIlnRLrIa8LlIaHx8HGfPnkVmZiZKS0upO097xjlsyHWsY/8vPyox+jh5ci7KDFNxuj9WAuMz8+jr68PS0pLNd7+7L8myuM2lfzljcZtf1JsegkfSMLgEocS4Px+5KE0acq4WAZJIJDh37hwiIyNRVlbm0FqIn58fEhISUFJSgs2bN6OiogIBAQEYGBjAsWPH0NjYiKGhIcp7zJ2IRCKcP38emZmZeg3gvs7MzAxaWlqowZe6xQ26a0lkcYOxKGm1puM8KhJyBNIKfnR0FBUVFXqlr+S0RlupTLM9hafLRYX6lWpkLwJ5LItyc4vPLKRFB0IsFmNoaAh+fn6Ii4tDfHy8VQPZflafice/MV+9NiNWYXBWjMxY0yf+7y/Jw6sNo2b3AwCbn2xA6wP6k0p9tQTbEqQNT2pqKu0XYhaLhcjISERGRiIvLw8ymQwzMzOYnZ2FQCBAUFAQ1ZPk6sF9QqEQzc3NyMvLQ1qa5TYAX2F2dpYSIN1IxtYScCYd5wHY+2NVqVTg8/mQyWSoq6tb8UXaO9juxWN9dh0PALzwI323bGMjGIQWzL3LshJRXp5FOTvPzMygtbUVWq0WsbGxiI+PR0xMjMkGz92XZuPA1+ZnFl3x/Dl0GIiHIe/trMANFtZ+NAC6JxdQkLgs3KQLdmZmJjIzM1fNHTFpw5Obm4v09HSnv15QUBDS09ORnp4OtVq94jyJiYmhRMmZaTvyQlxQUEANpFsNkO+7uLjYYirNXKMsuYadlWU5g+FreJQI2YNYLKby7uvXrzf6Q7NXhP7v6JBdx+QH4OI1SdR/W21AasCm/58O1HV2JghiRRUV6VkWFxenl+66uS4DT3zdD0vLovd+2I5HrjU9YqIkOQrpkYEYFplfI7v+pSa0PHARxsbG0NXVhTVr1iApKcnsc3yJyclJtLe3u+19c7ncFYP7ZmZmMDQ0hPb2dkRERFCCRKdR5vT0NFpbW1fd903eaBUVFSEx0Xx/niG6UZJCocDOnTsRHx+P22+/3RmH6tF4tQhNT0+jpaUF6enpyMvLM/mjsleE7O3U+OgXa6h/mxrBcLDJcoorLyF8xWOG6RjSs2xychLd3d0ICwtDfHw8NSH05Z+W4hdvtZo/3tZp7LosD7HBpht4v/xtHdbsP2J2P2oAvb29GBkZWVUl2AAoG57y8nLExsa6+3D0Bvfl5uZCLpdTabu+vj4EBARQguTI4D5SeEtLSxEfb/2YEm+HjHiLioocEl6VSoWbb74Zo6OjOHHihNud492BR4mQtXdmBEFgYGAAfX191EKgOewVIXu58tUOdO9L1jNKNIyAnjsyaHYf1t6jhoSEICQkBJmZmVAqlZiZmaEmhJIXmiAuILOgqJueOGkxLffA1lzs/7LX7DYn2kdw3Q9qVs0CK9n7NDIygsrKSkRGRrr7kIwSGBiItLQ0pKWl6aV3dQf3kaJk7ViF8fFxdHV1eYzwugqhUAg+n4/CwkKHBEitVuMXv/gF+vr6cPjw4VUpQICHiZA1aDQatLW1YX5+HrW1tVb1/7hahADoVb4YS8HNWpjKev1a20s1/f39qVEDpFXRzMwMHqoG/nDa8vN1I53CWH88cFUhKtK//2H8uDbNogjxFVG4eRUJUFdXF2ZmZlBdXe01wmuY3l1aWsLMzAxGRkbQ0dGB8PBwSpBMuTqTkV9FRcWqinjJ4ovCwkKHnM81Gg3uuOMOtLe348iRI6sqijTEq0SItDdhs9moq6uzegKhO0ToyUOduOvS5SY9e3LvD169xvJGZuBwONT6wJo1a/D28Fk0jUv//1/1+5GM0TWrxE9eb6H+O4QDnL33Qouvq1K79nN2F1qtFm1tbVhaWkJNTY1XOhoAy9mH8PBwhIeHIycnBwqFgkrb9ff3w9/fXy9tx+FwMDQ0hP7+fo+O/JzB/Pw8mpubUVBQ4LAA/eY3v8HZs2dx9OhRm9eTfA2vEaH5+XnweDzqompLDtteEarPikDDwILNzwOAt89O4veXFdn1XLphsVh4a+c6nUhHvx/JmuSfRAMUP3TU4nb12b6fliHnH6nVaj0bHl8gICAAqampSE1NhUajoQb3dXZ2QqVSUQ7gZWVlq06AeDwe8vPzHar+02q12LVrF44dO4ajR4+uqkpCU3hU44apiGFkZASNjY3Izc1FcXGxzYuo9orQazfX2P0ByTzQWeVftxiz4TEUJFtYuf31Vb49nI204QFAy1hqT4bD4SA2NhZFRUXYsGEDEhISIJfLERQUhObmZpw5cwb9/f1YXFz0CishexGJRJQApaZabt42hVarxR//+EccOnQI33zzzaqZpGsJj4uEdAfbabVadHV1YWJiApWVlXYv3DmSjuvcdwnu+GcjDveIbHpegIOfbPvYAopTHGuWNaQ0LQZliYFomTRVam1L2tB4BOVpFjJ0IpPJ0NTUhLCwMK9yhHYUgiAgEAgwNzeHdevWISQkhBrcNzMzg8HBQXC5XMTGxiIuLg7R0dFuGVjnDEgBysvLc1iA7r//fnz00Uc4cuQIcnJyaDxK78bjRIhEqVSiubkZKpUKdXV1Dk0cdHRN6IWbln3RCvZ+Y/VzsmPMrxEEsgG5mQaeHW824/w9m61+PWt595eWS62Noys61qXwfAmyHy0uLg6FhYWrpvmWIAh0dnZibm4O1dXV1O/QcHAfmbbr7u6GQqGgetdiY2O91qqJdL7Izc11yAGCIAj86U9/wjvvvIPDhw8jPz+fxqP0fjxShJaWltDU1ITw8HBa5tCTDWGO4scGVNYZ4iI9jIvjghmMz8shU6uREhmEDTkxCAlYjhS2Fsfjo9Zpk8+XqbSYW5IhJoz+Be/2+y+0an1HH10BMs6LN/jmj4u8G05PT18VhpwkWq0WHR0dWFhYMGs8y2azERMTg5iYGBAEQfWuTUxMoKurC6GhoZQghYeHe8XnR5og5+TkOCxAjz/+OF599VUcPnwYa9Y4VnDki7AID0vmDg8Po6WlBVlZWbT5bvX392NpaQnl5eUO7eemV8+hcdi+QgWS9FDgn79cj0A/LmoeP2F221A/Fs7ec6FDr2eOfzcOYu/n9k9HXWY5KuIC+O5/Kr3mImMtpC3LavND02q1aG1thUQiQVVVldWVqIYolUqqVWBubg5sNptK28XExHhk2m5xcRHnz59Hdna2Q+s2BEHg6aefxlNPPYVvv/0Wa9eupfEofQePEiGCIHD69GmkpKTQamk+NDSEubk5VFZWOrSfcZEUFz3d4PDxBHKA7/5wATY92QCZ2nxo9dGtlchPondtyBgEQeB41xQe/0qAwQXbvCJ+kB2KX5YHYnZ2Vq8HJTo62qvXTSYmJtDR0YHi4uJVVUar0WjQ0tIChUKByspK2oovtFotRCIRVQIul8upwX2xsbEeUeZOpwD99a9/xWOPPYZDhw6htraWxqP0LTxKhIDlOye6D2l0dBQTExOoqalxaD8EQaDwwW9pOaarSmKwoyYNN1gYtR3EAc7vMe9k4CxEIhEueY4HqYm/P35tHq4u/X6xVndtYHp6Gmq1mrrrdbaBJt0MDw+jt7cX5eXlq6qTXaPRoLm5GRqNBmvXrnXqd0am7WZnZyESiRASEkKdKxERES6PqJeWlqgxFJmZmXbvhyAIvPTSS3jwwQfxxRdfoL6+nr6D9EE8bk1ItzqOLhwtTCBHMGi1WsQGczErtddV7nuOC0Q48MMy+LNZUGpNv1+ZBuidWkSuER85Z0KmoQ7uWL4jtOaCoLs2UFBQQHXiDw4Oor293SvGVRMEgb6+PoyOjqKqqsquibzeilqtBo/HA4vFomUt1hK6llMqlYpK25HHoJu2c/axkAKUkZHhsAC9/vrr2Lt3Lz799FNGgKzA4yIhtVpNu7vB9PQ0BAIBNmzYYPNzdR2wAWBmSYKLnm2k5bg6HrgIHROLuO5l86O2E0L9cWSX7cduL6Ojo+ju7qY1DUXOvZmZmcH8/Dy1WB0XF4ewsDCPWEfSrQSrrKxcVbNdVCoVmpqa4Ofnh/Lycreu1Wi1WiwsLFAl4FKpVC9t50ilrDHIUeQZGRkOjVIgCAL//Oc/8fvf/x4ff/wxLr74YhqP0nfxuEjIGdgbCRmbAZQQEUbrsa1JCkeoPwdipenjmxIraX1NU5BRAGnGGRUVRdu+defeqFQqzM7OYnp6Wm9gX1xcnEOOzo5AehJKJJJVNYIc+L4BNygoyOEpsHTAZrMRFRWFqKgo5OXlQSqVUoLU09OD4OBgvbSdI8dLClB6errDAvTvf/8bu3btwsGDBxkBsgFGhEygK0CG/m9ZUQEYmHe85Jvk3V9U4qoXzpnd5uPmCWyrcN6sFq1Wi/b2dohEItTW1jo1CvDz80NSUhKSkpKg1WohFAoxPT2NtrY2amAfeZFxdhoGWI6+m5ubodVqUV1d7dMuCIbI5XI0NTUhNDTUYxtwg4OD9Qb3kWk7Pp8PAJQDeExMjE1rWGKxGOfPn0daWhqys7MdOsaPPvoIv/71r/Gvf/0LW7dudWhfqw2PEyFnpGVsFSFzAgQAF+bHYODMOG3H90Gz5X29e37EaSJETqbVaDSora21uxzXHsiS3djYWL1BbP39/WhrazM5sI8uFAoFeDwe/P39sXbtWo8sGXYWMpkM58+fR1RUFNasWeMRKVFLcLlcJCQkICEhgRrwODs7i4GBAbS1tSEyMpK6iTF3IyWRSHD+/HmkpKQ4LECffvopfvnLX+Ktt97CVVdd5dC+ViMeJ0LOwFoR0i1AMDcF9fZNGXjNQRGKDgCOd0/h7o+7IZJbPrZxEX2Rly6kM3lQUJDbL8KGg9ikUimmp6f1BvaRgmRqxIAtkBfhiIgIuzwJvRmpVIrz588jNjbWax0gdAc85ubmQiaTUWm73t5eBAUFUYIUGRlJfb8SiQSNjY1ISUlxuBfxyy+/xC233IK///3vuPbaa+l6a6uKVSNCpLiYutAYFiCYG8MdFuT4HblQAdz+rw6rt1+UO16Rt2Kfi4uUM3lBQYHHXYSDg4OpclnSq2x6ehoDAwPUwD7DC4y1kK4cCQkJKCgo8MqLsL2QaaikpCSzE4m9jaCgIGpwn1qtpgb3tba2QqvVIiYmBuHh4RgcHKRFgA4fPowdO3bgb3/7G2644QYa38nqwuNEyFnpOGB58dnYxcpYAYIlLsiKwAk7xzzYQyCX3giF/HGSTXmefiHS9SojJ4NOT0+jtbUVBEHYVM5LzoUhq6E8/b3TCVmKnJqaSpsjiSfC5XKpeVpkmndiYgK9vb0gCAIikQiDg4NU2s7Wz+H48eO48cYb8dxzz+Gmm27y2c/RFXicCDkDXREyXLi0Jv1mjJf+u9JOI1D7KE+lryrPGSXYrsRwMujCwgKmp6fR29u7Yh3JcH2LFF9Hbfm9EdIPLTMz06FKMG+DxWLBz88P09PTSEtLQ3p6OlXc0N/fj4CAAOomxprqzIaGBtxwww04cOAAbrnlFkaAHMTj+oS0Wi1UKvOjr+3hq6++woYNG/QWK8kIyFYBIlGr1Sh75Du6D9UoZ3dvQGigY1VbBEGgt7cXY2NjKC8vp7UE21OQSCSYnp7GzMwMFhcXqVHV8fHxWFhYQFdXF4qLi2m1hfIGyOgvJycH6enp7j4clyKTydDY2Ij4+Hjk5+fr/c7JqJp0blCr1VS1XWxs7IpKyXPnzuGaa67BQw89hN/85jeMANGAx0VCzvpSDYsTyOjHXgEC4JLyYQDwY8FhASJLsElHZF9txAwJCUFWVhaysrKoUdUzMzPo6+sDQRBISEhAQEAAVfm4GpibmwOfz1+V0R8pQHFxcSsECFgZVS8tLWF2dhYjIyPo6OhAeHg4BAIBcnNzweVysX37dtx///2MANGIx4mQsyBFyJYCBE/h67scs/7QLcGuqalxaQm2OwkICEBKSgpkMhlEIhGysrIgFovR3NwMANTFx1PdnOlgZmYGLS0tWLNmDZKSnNdn5onoCpA1xScsFgvh4eEIDw9HdnY2FAoFZmdn8fzzz+POO++ESqVCTU0N1qxZA4VCsaoamp3JqhIh0hLIlgIES1xaGIuvu2Yd3o8pLsgKR3yo/aLhSSXYrkar1aKzsxNCoVCvAVfXzbmnpwcKhYJKwcTFxflMs+rU1BTa2tpQUlKy6tKPcrmcKkG3t/qRvInZvXs3Dh8+jC1btiAqKgp33HEHZmdn8cYbb+C6665zwtGvLjxuTQgALQPoDGloaEBWVhZiYmJojX7kShUqHzM/F8heAllA0/32O2iTJdirsQxZo9GgtbUVMpkMa9euNXnXSg5hI9eRlpaWEBERQa0j0e1T5iomJibQ2dmJ0tJSxMXFuftwXIpcLkdjYyOio6NRVFTk0Hnf3d2Nyy+/HLfccgseeeQRymC5ra0NsbGxqy66dAarQoQIgkBjYyMUCgWSkpKQkJBAq4tzx/girnvFvAmprXAAtD5gvwB5Wwk2nahUKirlVlFRYZOVi1wup9aRhEIh5VMWHx/vNQP7RkdH0dPTs+rGUAD0ClBfXx+2bt2KH//4x/jzn//scX10voJHihCdM4XI9R+5XE41O5IuzmQfQWhoqMOvo1SpsPbRE2aGX1tPIBdoutd+ARoZGYFAIFiVVWAKhQJNTU0IDAxEWVmZQ+lHtVpNdeDPzs6CzWbrDezzxNTm8PAw+vr6UFFR4ZPVj+YgU3CRkZEO2xANDg7i8ssvx9VXX43nnnuOESAn4rMiRFrwkBVxuh5wKpWKGrw2NzeHoKAgSpAcHSvAG5rBTW+02f38/Bg/fPTrC+x67moowTaHVCpFU1MTdRGi88KhO7BvZmYGKpUKMTExiI+P95iBfQMDAxgcHERlZeWqmoMELN98NDY2UhZMjvyGR0dHcdlll2HLli144YUXGAFyMj4pQoYOCMZMSElIV15yTcDPz48SpMjISIfTL1qtFkKJAuFBfvjBk99hzkym8YUfFWJzgX05Zq1Wi7a2NiwuLmLt2rU+W4JtCtKGxxVWNARBQCwWU+eMWCxGZGQk4uPj3TKwz3AQX1gYveNGPB2FQoHz588jPDzcYQGamJjA1q1bsXHjRrz88sseGe36Gh4pQiqViiqhthWyBJu06LHlhNRqtXqCxGKxqPWA6OhoWu6Ibnr1NHhjMr3H8mMD8NGv7C/DJtdAtFot1q5d6zPVXdYiFArB5/ORlZXllvUvw4F95JhqOiJrSxAEgZ6eHkxOTqKqqoqW1LI3oVQq0djYiLCwMJSUlDj0WU9NTeHyyy9HdXU13njjDUaAXIRPiZCjDgi6kGW809PTmJ6ehkajoS4sntRXQpZgBwcHo7S01GOOy1WQc4gKCgqQkpLi7sOhBvaR60hcLpc6b+ge2EcQBLq6ujA7O4uqqiqvreSzF10BctQFfXZ2FldccQXWrFmDt99+22WN6Aw+JEJ0OCCYgjRAnJ6extTUFBQKBWJjY92+HrCwsIDm5uZVWYINAGNjY+ju7kZJSQni4+PdfTgrIAf2kVGSRqOh1pFsHcBmCEEQ6OjowPz8PKqqqlyeAnQ35DTYkJAQh4fxCYVCXHnllcjKysJ777236jIJ7sYjRYhsKrUGcwUIzoDsK5mamsL09DQkEgmio6OpdSRXncBkCTbpBbaaBIggCAwODmJwcBDl5eWIjo529yFZRHdg38zMDCQSCaKioqh1JFu678n1P7FYjMrKylXXuU8KEBn9OyJAIpEIV199NRITE/HBBx+sGjcRT8KrRciWAgRnQQ5em56exuLiIiIiIpCQkODUBerVXIJNEAQEAgEmJiZQWVnptYvwUqmUEiSRSITQ0FAqbWduYJ9Wq0VLSwtkMhmqqqpW3V27SqXC+fPnERQU5LAALS4uYvv27QgPD8cnn3yy6sTcU/BaEdJd/2GxWB5RRkk2OpK9SGFhYVSEREe1mm4JdkVFBSIjIx0/aC9Cq9Wio6MDIpEIlZWVPrMGQg7sm5mZwdzcHFWhaTiwT6PRgM/nQ6VSobKy0iPKwl0JKUBkD5gjv3mJRIIf/vCH4HK5+Oyzz3zmXPJGvFKE6CxAcBbkhWVqagpCodDhXiSNRoP29vZVW4Kt0WjQ0tICuVyOyspKn02b6I4WmJmZgVarpZpjx8bGAABr165ddQvndAqQVCrF9ddfD7VajS+++GLVVRR6Gh4pQhqNBmq18XHW3iBAhpCd99PT05idnaXudBMSEhAREWHxPZAl2ARBoKKiYlWmYOy14fFmyIF9k5OTGBsbo0ZUk1GSrwqxISqVCk1NTfD390d5eblDAiSXy/GjH/0IYrEYX3755apr6vVEvEaEXF2A4Cx0R1OTvUhkhGSshFcqlYLH4yE0NBQlJSWrrgRbLpdTLuCrsQRdqVSiqakJAQEByM3NpSaCLiwsUAP77B1R7Q2o1Wo0NTXBz8/PYQFSKBS46aabMD09ja+//nrVOYp4Kl4hQoYzgLxVgAyx1ItEzr5JTEw0OpDL15FIJGhqaqLMKD1h3c+VkE4A5A2I7vtXKpXU+qNQKERAQIDeOpIvnCukAHG5XJSXlzt0A6JSqbBjxw4MDQ3h22+/XXXGrp6MR4qQ7ohvwwo4X70QkakXUpAUCgW0Wi0SExNRVFS06tYAFhcX0dTUhJSUFOTm5vrERdUWSDPOiIgIiz54Go2GipBmZmYAeP/APrVaDR6PBw6H47AAqdVq/PznP0dXVxeOHDmy6kZbeDoeLULeuP5DB8PDwxAIBIiNjYVUKoVEItFbC/D1NSFdG57MzEx3H47LkUqlOH/+PGJiYmweR0AQBDWwj7yZ8baBfaQAsdlsVFRUOCRAGo0Gt912G3g8Ho4cOYLExEQaj5SBDjxWhJRK5aoTILIHZnx8XK8E27AXiTTLjI+P97neBnIaaFFREZKTk919OC5HIpHg/PnzSEhIcDgFSzZWk4KkO7CPXEfyNDQaDZqammgToN/85jc4efIkjh496hG2Tgwr8UgR+vTTT6FWq7Fx40YEBgauCgGytgSb7EWampqCSCSivRfJnZDD2FbjNFDgeyfwlJQU5OTk0H7e687U0h3YFxcXZ1WVprPRaDTg8XgA4PAoeq1Wi9/97nc4fPgwjhw5goyMDFqO8a233sJ3332H8+fPo7W1FUqlEq+99hpuvvnmFds++OCD2Ldvn9H9BAQEQC6XG/3b22+/jWeeeQbt7e3w9/dHXV0d9u/fj+rqalreg6fhkQsNXV1dePrppyGVSnHllVdi27Zt+MEPfuBzd/0kSqUSfD4fBEGgtrbWbMokMDAQaWlpSEtL01uc7uvrQ0hIiN6gPndfVKyFIAgMDAxgaGgIa9euXZVVSwsLC+DxeEhPT0d2drZTXiMwMBCpqalITU2lRpjMzMxQqS93DuzTaDRUG0JlZaXDArR792589dVXOHr0KG0CBAD33XcfhoaGqNHeQ0NDFp/zs5/9bEVa2dQa7yOPPII9e/YgPT0dt99+O8RiMd59911s2LABhw4dwoUXXkjDu/AsPDISApZPyoaGBrz//vv48MMPIRKJsHXrVmzfvh2XXnqpz3Q401WCbdiL5O/vTwmSJ9zlmoIgCHR3d2NqasqrbXgcQSQSgcfjUaPYXQ1ZpUne0JAD++Li4hAbG+v0dSRSgMhRJI4U4Wi1Wtx3333497//jaNHjyIvL4/GIwW++eYb5OXlISMjA4899hjuuecei5HQkSNHrBIPgUCANWvWIDs7G2fPnqV6mNrb21FbW4ukpCR0dXX5XJGSx74bDoeDjRs3YuPGjXjqqadw7tw5HDx4EPfffz9uvfVWbNmyBdu2bcPWrVu99sJF3v0mJSU5nP/ncrlITExEYmKiXi8SeZdrrhfJXWi1WrS3t2NhYQG1tbWrzgkaWC7CaG5uRl5eHtLS0txyDGw2G9HR0YiOjkZ+fj41sG94eBgdHR2IjIykoiS6b/5IKyI6BIggCDz00EN49913nSJAAHDJJZfQvk+S1157DWq1Gnv27NFroi0uLsaOHTvw4osv4vDhw9iyZYvTjsEdeKwI6cJms7Fu3TqsW7cOjz/+OJqbm3Hw4EE8+uijuP3223HJJZdg27ZtuOKKKzz6rl8Xcg5Obm4u0tPTad03h8OhLhpFRUWYn5+nXo8gCL1Bfe4q3yUvPkql0mIK0leZnZ1FS0sLCgsLPaYIg8ViISwsDGFhYcjJydHzQxQIBNTAvri4OISHhzv0WyPPAbVajcrKSocF6LHHHsPf//53HD58GIWFhXbvi26+++47nD17FhwOB4WFhbjkkkuMul0cPXoUAIyKzGWXXYYXX3wRx44dY0TI3bDZbFRWVqKyshIPP/ww2tracPDgQTz33HP49a9/jYsuugjbt2/HlVdeiejoaI8UpOHhYfT29rpkDg6bzUZMTAxiYmJQWFhI9SJ1d3dDqVTqzUVyVZivVCrR3NwMNpuN6upqn0svWMP09DRaW1tRXFzs0WXDumuQKpWKmjzc1NSkd7Nj6+Rh0g2cLgF6+umn8fzzz+Pbb79FaWmp3ftyBg888IDefyclJeGNN97ApZdeqve4QCBAaGio0fOBjOoEAoHzDtRNeEZexk5YLBZKS0uxb98+8Pl88Pl8ajZ8Tk4OrrnmGrz66quYnp6GJyx9kaOY+/v7UVlZ6fJBbCwWC5GRkcjPz8eGDRtQU1ODkJAQDAwM4NixY+DxeBgbG4NSqXTaMcjlcjQ2NiIgIGBVGnECwMTEBFpbW1FaWurRAmSIn58fEhMTUVZWhs2bN6O4uBgA0NHRgaNHj6KlpQUTExNUo7kptFotFQXTkYL761//iieffBJffvklKioq7N4X3VRUVOCNN97A4OAgZDIZBAIBHnroIYhEIlxzzTXg8/l62y8sLJj0sgsPD6e28TU8tjDBEQiCQH9/P95//3188MEHOH/+POrq6rB9+3Zcc801SEpKcnmEpNFoqEFka9eu9bjCColEQvUiLS0t2T1wzdJrNDU12dWE6SuQ02DLy8t9xjqGIAgsLS1RfojkwD4yStJd6yMjIIVC4fA4CoIg8NJLL2Hfvn344osvUFdXR8fbsRpLhQmmePnll/HLX/4S1113Hf79739Tj5PFRKOjoyueMzIygvT0dGzZsgWHDh2i4/A9Bp8UIV0IgsDw8DA++OADfPDBBzh9+jRqampwzTXXYPv27UhLS3P6xZBMPwHwChdsuVxOCZJIJEJ4eDhV2GCveJJFGKmpqU7pgfEGyDRsRUWFV0yDtReZTEatI+kO7IuNjUV/fz8UCgWqqqocFqDXX38d99xzDz799FNs2rSJxndgHfaKkFKppNbWxsfHqcfj4uIgl8uxtLS04jnt7e0oKSnB9ddfj/fee4+Ow/cYfD4XwmKxkJGRgbvuugu/+93vMD4+jg8//BDvv/8+7r//fpSXl2P79u3Ytm0bsrOzab84eqMLdmBgINLT05Genq7Xi9Tb22tXL9Lc3Bz4fL5TijC8hcHBQQwMDKCystLnhxEGBQVR549KpaLOn4GBAQBAcnIyFhcX7a7UJAgCb731Fu6++2588sknbhEgR/D390dYWBikUqne43l5eTh16hQmJydXpGnJtSBnVPy5G69eE7IVFouFlJQU3HnnnTh8+DBGR0dx66234vjx46iqqkJ9fT0ef/xxdHV10bKGtLCwgLNnzyI2NhZlZWVeIUCG+Pv7IyUlBWvXrsXmzZuRmZkJsViMs2fP4uTJkxAIBFhYWDD5eU1OToLP56OoqGhVChBBEOjr68Pg4CCqqqp8XoAMIdeRWCwWgoODqXWktrY2HDt2DK2trZicnDQ5P8wQgiDw3nvv4fe//z0OHjyIiy66yJmH7xQEAgHm5+dXNLBu3rwZAPDVV1+teA6ZgiO38SV8Ph1nDQRBQCgU4uOPP8YHH3yAb775Bjk5Odi2bRuuvfZau8YIOLME2xMgnZvJdQAOh0NFSORI6pGREQgEApSVlSE2Ntbdh+xyyHHs4+PjqKqqWpUTPLVaLdra2iCRSFBVVUWlogmCwOLiInX+SKVSREdHU+tIptYhP/jgA9x2223417/+hauuusqVb2UF5tJxS0tLGBgYQFlZmd7j8/Pz2LZtG7777js89thj+OMf/0j9raenB8XFxauuWZURISOIRCL85z//wQcffIBDhw4hNTUV27Ztw/bt260arOXKEmxPQKvVUr1IZCViYGAgpFKpz69/mIJ0gpiZmUFlZaXX+/rZA9mMLBaL9QTIGKTRKjmwLywsjBpFERYWBjabjU8//RS33HIL3nrrLVx77bUufCff88orr+DEiRMAgNbWVjQ1NWHDhg3Izc0FAGzfvh3bt2/H4OAgsrKyUF1djdLSUsTHx2NsbAxffPEF5ubmcOmll+LTTz9d8Zk8/PDDuO+++5Ceno7rrrsOEokE77zzDmQyGQ4dOuSVkZ8lGBGywNLSEj7//HO8//77+OKLLxAbG4trrrkG1157Laqrq/UEiSzBnpiYwNq1a1fl6GCy+kkoFILL5UKtViM2NhYJCQmIiYnxubs4YxAEgY6ODszPz6OqqmpVOkEQBIG2tjYsLS2hurrapmIcch1yZmYG77//Pt555x0UFhaioaEBr732Gm688UYnHrl5br75Zrzxxhsm/7537148+OCDWFxcxL333ovTp09jaGgIIpEIISEhKC0txU9/+lPs3LnTZHr+n//8p1ED05qaGme9LbfCiJANSKVSfPnll3j//ffx2WefITw8HFdffTW2b9+OkpIS3Hzzzbjyyivxk5/8xONKsF0BmXpZWlpCZWUlAgMDqdLd6elpyGQyvblIjlRHeSrk3f/i4iKqqqp81nTXHARB6H0GxtwBrGVpaQmPP/44nnvuOYSEhMDPzw9XXnklrr32Wmzfvp2+g2ZwG4wI2YlcLsfXX3+N999/Hx999BEkEgnCw8Px7LPP4pprrlkVd/y6qNVqyoJl7dq1Ru98TfUixcfHO3Sh8hS0Wi1aW1shlUpRWVnpE+/JVkgBWlhYQHV1tcOfwfHjx3H99dfjueeew3//93/jzJkz+PjjjzE8PIx3332XpqNmcCeMCDlIb28vtm7diqSkJOTk5ODzzz8HQRC46qqrsH37dmzevNnj+4IcRalUgsfjgcvlory83CoBJntJpqamsLCwQEsvkjshfdBUKpXDTZjeCpmGFIlEtESBDQ0N+OEPf4gnnngCt95666rsLVsNMCLkAAqFAvn5+bjhhhvw+OOPg81mQ61W4/jx4zh48CA++ugjyGQyXHXVVdi2bRsuvvhin0vPyGQyNDU1ISwsDCUlJXb1fSgUCqqXRCgUIjQ0VG9Qn6dffNRqNTULp6KiYtUKUGdnJ4RCIaqrqx0+z8+ePYvt27fjoYcewp133unx5wCD/TAi5CD9/f0mh5BpNBqcPHmSmom0sLDgUzORxGIxmpqaEBcXh8LCQlouFCqVSm8uUmBgICVIjro2OwOVSgUejwcOh+PwOGpvhW4BampqwtVXX4377rsPu3bt8rjvnIFeGBFyEVqtFmfPnsXBgwfx4YcfYmpqClu2bMH27duxdetWr+shEYlEaG5uRlpamlOcJgDrepHciVKpRFNTEwICAry2GdlRCIJAV1cX5ubmaBGglpYWXHnllfjDH/6Au+++mxGgVQAjQm5Aq9WCx+Ph4MGD+OCDDzAyMoIf/OAH2L59O6644gqPvOPXhZyD48pBbFqtlhrUNzMzozcXKSYmxuWCpFAo0NTUhODgYJSWlrpdEN2Bbi9UdXW1w6XoHR0duPzyy3HnnXfigQce8OjfAAN9MCLkZsh+in//+9/48MMP0dPTg4svvhjbtm3DVVddhaioKI/6MU5MTKCjo8Otc3AIgoBIJKIq7VQqFSVIsbGxTo9I5HI5zp8/j/DwcBQXFzMCRIMAdXd34/LLL8fPf/5zPPzwwx51zjM4F0aEPAgytUGOoGhra8OmTZuwbds2XH311YiLi3Prj5N0gvCkMQS6YwSmpqYgl8ud2oskk8lw/vx5REdHr9pxFGRT9vT0NC0C1Nvbi8svvxw33ngjDhw4sCpFfTXDiJCHQhpfkoLU1NSE+vp6bNu2zeUzkchjGR0d9XgnCLFYTEVIYrGY8iOjoxdJIpHg/PnziI+PR0FBwaoWoKmpKVRXVztcXDM4OIitW7di27ZtePbZZxkBWoUwIuQFkDORSEE6ffo0amtrsW3bNmzbts2pM5HIyqe5uTmsXbvWqwooZDIZJUjk1EqysMHWu3exWIzz588jOTkZubm5q1aABAIBJicnaRGgkZERXHbZZdi6dSuef/55RoBWKYwIeRkEQWB8fBwffPAB3n//fZw8eRIVFRXUTKSsrCzaLpCkA4BEIqFseLwVc71IloR1cXERTU1NTq0E9HR0HcGrq6sdNmSdmJjAZZddhk2bNuHll19elZWFDMswIuTFEASBqakpfPTRR3j//fdx7NgxrFmzhnL8zs/Pt/uCSdrwaDQar5gGawu6g9bm5ubM9iKJRCLweDxkZWWtmP+yWiDTsWNjY7QI0NTUFC6//HLU1NTg9ddfZwRolcOIkI9AzkT66KOPqJlIeXl5lOO3LTORyP4Xf39/lJWV+bQPnkaj0WuO5XK5lCBptVrw+XyXlqJ7Ir29vbQJ0OzsLK644goUFxfjn//8p0+fWwzWwYiQD0IQBBYWFvDJJ5/ggw8+wFdffYW0tDRKkMrKykwKEmnDsxrLj3V7kaampqBWqxEZGYmsrCxER0evqs+ChCxIoWMon1AoxJVXXons7Gy89957q9LeiGEljAitApaWlvDZZ59RM5Hi4+MpQaqqqqIurr29vRgdHUViYuKqrf4ClqfitrS0IDMzE2q1GtPT09BoNIiNjXVZL5In0N/fj+HhYVRXVzssQCKRCFdffTWSkpLwwQcf+FR6l8ExGBFaZUgkEr2ZSBEREbjmmmuQnJyMRx55BC+++CKuvfbaVStAk5OTaG9vp6ZhAvqjqKenpyGXy/UEyRfv6AcGBjA0NISqqiqEhYU5tK/FxUVs374dERER+Pjjj726wIWBfhgRWsXIZDJ8/fXXePbZZ3H48GGEhYXhxhtvxPbt27Fhw4ZVl68fHx9HV1cXysrKEBsba3QbgiD05iKRvUhkc6wvzBCiU4DEYjF++MMfwt/fH59++qnXm/Yy0M/qusow6BEUFASRSITTp0/jnXfeQWRkJA4ePIgdO3aAxWJREyw3bdrk8+mTkZERCAQCVFRUIDo62uR2LBYLoaGhCA0NRXZ2NqRSKaanpykBc6QXyRMYHBykTYCkUiluuOEGsNlsfPLJJ4wAMRhl9a20MlAoFAo899xz+OSTT/DjH/8YW7duxSuvvIKJiQm88847CAwMxG233Ybs7Gzcdttt+OKLLyCXy9192LQzNDSE3t5eVFZWmhUgYwQHByMzMxO1tbW44IILkJiYiNnZWZw8eRKnT5/GwMAAJBKJk46cXgYHBzEwMIDKykqHBUgul+PGG2+EUqnEp59+SluT81tvvYXbbruNmtrKYrHw+uuvm9x+cXERu3b9v/buPS7mdI8D+GeSiu5ohW50s3LpoojIPV1O5cguh80ti2VpYsPmpOyWtEfscXblfg0vr1gskkgoJZcMDSociq4uU0qlmuf8sWfmZRSimaap7/v18nrl93v69TzT5TvP8/v+nm8AjI2NoaqqCmNjYwQEBKCsrOy9n7N//344ODhAXV0durq6cHNzw7Vr16TSf1IfLce1cUKh8INZX6KaSKISFGVlZXB1dYW3tzfGjBmj8O9uRTffbW1toaWlJbXrvvssUocOHcQzJE1NzRZ3z+3x48d4+PAh7Ozsmvw6VFdXY+rUqSgpKUFCQgJ0dHSk00kAJiYmePz4Mbp06QJ1dXU8fvwYO3fuxIwZM+q1raiogJOTE27evImxY8fC1tYWPB4Pp0+fhrW1NZKTk+ulnIeHhyMoKAhGRkbw8fFBeXk5Dh48iKqqKsTHx2PEiBFSGwv5CwUh0mhCoRBXrlwRB6Ti4mK4uLjA29sbLi4uCrWlz9s7AEjjnf+H1NbW4vnz5ygqKsKzZ8/Qvn17ibpI8g5Iubm5ePDgAWxtbZu8L+CbN2/g6+uLvLw8nD17Vuob3YqefzM2NkZERARWrFjx3iC0atUqrF69GoGBgVi7dm2948HBwQgNDRUfz8nJQZ8+fdCrVy+kp6eLXws+nw8HBwd069YN9+7da3P3SmWNghD5LEKhEDdu3BDXRHry5AnGjBkDLy+vFl8TSVSGoLi4GHZ2dk1+APNT1NXVSdRF4nA40NPTQ9euXaGrq9vszyLl5eWJlyKbGoBqa2sxa9YsZGVlITExEXp6elLqZcM+FIQYYzAwMEBZWRkKCwslvsdVVVXo3r07OnbsiLy8PPHP6Y8//og1a9Zg9+7d8PX1lbje/PnzER0djfj4eIwbN06m42pr6J4Q+SxKSkoYOHAgIiIicO/ePaSlpWHAgAFYv349TExMMGnSJOzZswcvXrxAS3qfI9qQtaSkBPb29s0agACgXbt20NPTg5WVFYYPHy4uiMfn83HhwgVkZmaKn0uSNVEAksbO6LW1tZg7dy7u3LmDhIQEmQegj8nJyUF+fj6GDh1a73uspqaG4cOH4+nTp7h//774eFJSEgA0GGRcXFwAABcuXJBdp9soCkKkyZSUlNC/f3/89NNPuH37NjIyMjB48GBs3rwZvXr1gpeXF3bs2CGuiCovQqEQfD4fL1++hL29vdyz15SUlNCpUyf07t0bw4YNg62tLVRVVZGdnY2kpCTweDwUFBSgpqZG6l/7yZMnyMnJgY2NTZPv2dTV1eH777/H9evXcfbsWbkVO3xbTk4OAMDc3LzB86LjonaijzU0NBrsf0PtiXS06iAkEAiwaNEiODo6Ql9fH6qqqujRowdGjRqFw4cPN/gHkbJpmobD4aBPnz4IDg7GjRs3wOfzMWbMGOzZswdmZmZwc3PD5s2bUVBQ0KwBSbQj+KtXrzBw4MAW98Akh8OBtrY2zM3NMXToUDg4OEBDQwOPHj3ChQsXcOPGDTx58gRv3rxp8td68uQJsrOzpRKAhEIhuFwukpOTcfbsWXTv3r3J/ZOG0tJSAHjvDE+UfCFqJ/r4U9oT6WjVQejZs2fYsWMH1NXV4e3tjSVLlsDV1RV8Ph8+Pj6YO3euRPuKigo4Oztj/fr1sLS0BJfLRZ8+fbB+/Xo4Ozs3mGobHh6OqVOnoqioCPPmzcNXX32FlJQUDB06VDy9b6s4HA7Mzc2xfPlyXLlyBTk5OfD09ERsbCwsLS0xbtw4/Oc//0FeXp5MA1JdXR14PB4qKythZ2fX4h8o5XA40NTUhKmpKRwdHTFkyBB06tQJ+fn5uHjxIq5du4bc3NzPSpd/+vSpOADp6uo2qZ9CoRCBgYFISEjA2bNnYWRk1KTrkbapVad59OzZEwKBoF42y6tXrzB48GBs3boVixcvhpWVFQAgMjISN2/efG82TWRkZL1smlWrVsHCwkIim2bRokVwcHCAn58fZdP8H4fDgYmJCZYsWYKAgAA8ffoUR44cwZEjRxAUFAQbGxtxkT5p1kSqq6vDzZs3UVdXBzs7O4XcYkf0LJKJiQmqqqrEqd/Z2dnQ1NQUZ9p97P5Wfn4+srKyYG1tLZUAFBQUhOPHj+P8+fPo2bNnk64nbaLfxffNXEQrG2/PfLS1tT+pPZGOVj0TateuXYMBQFNTU3yjUXRjkjGGbdu2QUNDA8HBwRLtV6xYAV1dXWzfvl3iHfvOnTtRW1uLoKAgiR9OKysr+Pr64sGDB0hMTJTF0BQah8OBgYEBFi1ahPPnzyMvLw8zZ85EUlISbGxs4OTkhMjISGRlZTVphlRbW4sbN26AMQZbW1uFDEDvUlNTg6GhIezs7DB8+HAYGBhAIBAgNTUVly9fxv379/Hq1at6r5toR4eP7QjRGIwxrF69GocOHRKnTLc0H7uH09A9I3Nzc5SXl6OwsLBR7Yl0tOog9D5VVVVITEwU378AKJtGXjgcDvT19TF//nycOXMGBQUFWLhwIdLT0zF48GAMGjQIYWFhuHPnzicFpJqaGly/fh3t2rWDjY1Nq5yNqqiooEePHrCxscGIESPE2whdvXoVKSkpyM7OhkAgEAegAQMGSCUArVmzBrt27UJCQgJ69+4tpdFIl7m5Obp3746UlJR6y+hVVVW4ePGiuFS7iLOzMwDgzJkz9a4XHx8v0YZIT5sIQgKBACEhIQgODsa8efNgYWEBHo+H4ODgeu+YKJtGfjgcDrp06YLZs2fj5MmTKCwsRGBgIDIzM8XZYyEhIeDxeBAKhe+9zps3b3Dt2jWoqqrC2tq6TZRdUFZWhr6+Pvr37w9nZ2dYWlqipqZGnBwiWn770Ov2MYwxREVFYdOmTUhISEDfvn2l1X2p43A48PPzQ3l5OVavXi1xbs2aNXj58iX8/Pwkln1nzpwJZWVlhIWFSSzL8fl87NmzB6amphg1alSzjaGtaH1vDxsgEAgk7uW0b98ev/zyC5YsWSI+9rnZNKLt/hvTnjQeh8OBrq4ufH194evri7KyMnFNpLFjx+KLL74QlzF/uyZSfn4+Hj58CC0tLfTt27dNFqITPYtUV1eHgoICmJmZoaqqCpmZmRAKhdDT08MXX3yBzp07NzpAM8awceNGrF+/HvHx8RgwYICMR9Gwbdu2ITk5GQBw+/Zt8THRqoS3tze8vb0BAIGBgTh+/DgiIyORkZEBOzs78Hg8xMXFwdraGoGBgRLXtrCwQEhICFauXIn+/fvDx8cHFRUVOHDgAGpqarB169ZWOaOWtzbxipqYmIAxhrq6OuTl5eHgwYMICgrC5cuXcejQIfrBUgBaWlqYMmUKpkyZgoqKCsTFxeHw4cPw9PSEjo4OPD09YWtri6CgIMyfPx8BAQEtdseG5lBUVAQ+nw9ra2txWYrevXujtLRUnNTw5s0bibpI7/s9YIxh8+bNiIiIQFxcHOzt7ZtzKBKSk5Oxe/duiWMpKSlISUkB8NfvuigIqaurIykpCaGhoYiNjUVSUhL09fXB5XKxatWqBhM5goKCYGJigg0bNmDTpk1QUVHBkCFDsHr1armOuzVrU39927VrBxMTEyxfvhzt2rVDYGAgtm7divnz51M2jQJRV1eHj48PfHx8UFlZiTNnzmDnzp347bffoKqqitzcXFy6dAlDhgxpk28wioqKkJmZWa8uEofDgY6ODnR0dMQ34YuLi/Hf//4XmZmZ6Ny5s7gukqh0B2MMO3fuREhICE6ePAlHR0d5DQsAsGvXrg/umv0ubW1tREVFISoqqtGfM3XqVEydOvUzekc+R9tbq/g/UTKBaBpP2TSKqUOHDjA3N0d6ejoWL16M2NhYCIVCfPPNNzAzM8PChQtx7tw5qTzkqQiKi4vFAehDW+e8+yySo6MjdHR08OTJE1y8eBHffvstwsPDERUVhRUrVuDYsWMYNmxYM46EtBVtNgjl5+cDgPidMmXTKK7AwEDMmTMHUVFRcHd3x/bt25Gfn4/9+/dDRUUFc+bMQa9evTBv3jzExcWhurpa3l2WieLiYty+fRv9+vX75L3b1NXV0bNnTwwaNAhOTk4wMzNDTEwMgoODYWBggKtXr0pkhhIiLa06CN28ebPB5bIXL17gxx9/BAC4uroCoGwaRRYbG4vQ0FCJ70379u0xZswYREdH4+nTp/jjjz+gra0Nf39/9OzZE7Nnz8bx48fx+vVrOfZcekpKSsQB6H3JMo2lpqYGMzMzlJSUYN++fVi6dCkuXLgAKysrRERESKnHhPylVZdy8Pf3x7Zt2zBy5EgYGxuLi2CdPHkS5eXlmDhxIg4dOiTOoHq3CNa72TQNFcEKCwvDypUrxUWwRNk0lZWViI+Px8iRI+UxdPIeQqEQaWlpiI2NxdGjR1FSUgIXFxd4eXkpXE0kkZKSEty6dQt9+/ZF165dm3y9P//8E7NmzUJMTIz4Jj/w1/3S6urqJgc5Qt7WqoNQcnIytm/fjrS0NOTn5+P169fo1KkTbG1t4evri8mTJ9fLoCotLRVn0xQWFkJfXx8+Pj5YtWrVe5MMYmJisGHDBvD5fKioqMDR0ZGyaRSAUCjE9evXxUX6njx5grFjx8LLywuurq4KkVTy7Nkz8Hg8qQWguLg4TJ8+HTt37sSkSZOk0ENCPqxVByFCGksoFOLWrVs4fPgwjhw5ggcPHmDUqFHw8vKCu7s7dHV1W1zKtygAWVlZSaV8wrlz5zBlyhRs2bIFU6ZMaXHjJa0TBSFC3iEqfCeqGnvnzh04OzvD29sbHh4e6NKli9z/QD9//hw8Hg9ffvklunXr1uTrXbx4EZMmTcLGjRsxffp0uY+PtB2tOjGhNYqMjASHwwGHw0FaWlqDbagmUtO8XRMpIyMDfD4fo0aNwq5du2BmZgZ3d3e51EQSkXYASklJwVdffYV169ZRACLNjmZCCuTu3bvizTgrKiqQmpqKwYMHS7R5N7nC1tYWPB4Pp0+ffm9yRXh4OIKCgsTJFeXl5Th48CCqqqoQHx+PESNGNOMoWy7GGB49eiReshNtsurp6QkvLy8YGBjI/A/4ixcvcPPmTfTu3VsqBeTS09Ph5eWFsLAwLFiwgAIQaX6MKITa2lpmb2/PHBwc2LRp0xgAlpqaWq9dcHAwA8ACAwMbPB4cHCxxPDs7mykrKzMLCwsmEAjExzMzM1nHjh2Zqakpq6mpkc2gFJhQKGS5ublsw4YNbPjw4axdu3bM3t6ehYWFsczMTFZeXs4qKiqk+i8vL4/9+eefLCcnRyrXu3TpEtPR0WFRUVFMKBTK+yUlbRQtxymItWvXgsfjYceOHe/ddJJRTaRmw+FwYGhoiMWLFyMpKQl5eXmYMWMGEhMTYW1tDScnJ/zyyy/Izs6WypLdy5cvcfPmTVhaWkplBnTr1i14enpi2bJl8Pf3pxkQkRsKQgogMzMToaGhWLlypbgKbEOoJpJ8cDgcdOvWDd999x0SEhKQn5+PBQsWIC0tDYMGDcLgwYMRHh7+yTWRRF6+fImMjAxYWFigR48eTe7vnTt34OHhAX9/fyxbtowCEJErCkItXG1tLWbMmIEvv/wSy5cv/2BbqokkfxwOB3p6evDz88OpU6dQWFiIpUuX4tatWxg2bBjs7OwQGhr60ZpIIgKBQByADAwMmty/rKwseHh4YO7cufjnP/9JAYjIXdvbYljBhIeHg8fj4cqVKx8tT001kVoWUU2k6dOnY/r06SgrK8OJEydw+PBhjBkzBvr6+vD09MSECRNga2tbr/aRKACZm5tLJQDdv38fHh4e+Oabb/DTTz9RACItAs2EWjAej4eff/4ZS5cuha2trby7Q5pIS0sL//jHP3D48GEUFRUhIiICBQUF8PDwgJWVFZYtW4bU1FTU1dUhKSkJP/zwA8zMzGBoaNjkr/3o0SN4eHhg4sSJWLt2bZss9kdaJvpJbMGmT58OU1NThISENKo91URSHBoaGpg0aRIOHjyIoqIi/PrrrygtLYWPjw969eoFb29vKCkpSeU5oLy8PLi5ucHNzQ0bNmygAERaFPppbMF4PB7u3bsHNTU18QOqHA5HXFnS0dERHA4HR48eBUA1kRRVhw4d4O3tjT179uDkyZOorKzEgAEDEBcXBzMzM3z//fc4d+4campqPvnaBQUFcHd3x+jRo/Hbb79RACItDt0TasFmz57d4PGLFy8iJycHnp6e0NPTg4mJCYD6NZHezpD7UE2k1NRUnDlzBr6+vhJfh2oiNS8ejwcPDw+EhoZiyZIlqKmpQVJSEg4fPow5c+agpqYGHh4e8PLywsiRI6GqqvrB6xUWFsLNzQ1DhgzBli1b3pvaT4hcyfcxJfI5pk+fLrWHVbOysuhh1RaiqKiIbd++vcFzNTU17Pz582zBggWsR48eTFtbm02ePJkdPHiQPXv2rN6DqI8ePWJ9+vRhX3/9NX3/FNjEiRMZgPf+XDDG2MqVKxkAtnTp0mbsmfRQEFJAHwpC5eXlzNramgFgY8eOZcuXL2eurq4MALO2tmbl5eX1Pufnn39mAJiRkRELCAhgc+fOZVpaWqx9+/YsMTGxOYZEPkFtbS1LTk5m/v7+zMTEhGloaLCJEyeyvXv3sqKiIpabm8v69evHJkyYwN68eSPv7pImKCkpYV27dmVaWlrs8ePH9c5fu3aNKSsrsz59+rCqqio59LDpKAgpoA8FIcYYEwgEjMvlMkNDQ9a+fXtmaGjIuFyuxEznXfv27WMDBw5kHTp0YNra2mz8+PEsPT1dVkMgUlJXV8euXLnCfvjhB2ZmZsbU1NSYjo4OGzduHKuurpZ394gUHDt2jAFgo0ePltheqaqqillZWTFlZWV2/fp1OfawaSgIEdJKiALSuHHjGpzxyoKxsTED0OC/uXPn1mtfWlrKuFwuMzIyYioqKszIyIhxuVxWWlraLP1VVDNmzGAA2MaNG8XHli1bxgCw0NBQOfas6WgXbULIZzMxMYFAIIC/v3+9cwMHDoSHh4f4/5+zwzv5S1lZGfr16ycuZPj8+XMMHToUNjY2SE1NhbKyAueYyTsKEkIUl7GxMTM2Nm5U209NmiGSEhMTGYfDYY6OjszS0pKpqamxu3fvyrtbTUZBiMhccyzZxMTEMHt7e9axY0emo6PDXF1d2dWrV2U5LMIaH4SEQiHr3r0709DQqLdUWFlZyXR1dVmPHj2opMRHLFq0SPy7s379enl3RypoOY7InKyXbKgon/yYmJiguroaERERePr0KXR1dTFkyBAMGDBAol12djYsLS3h4uKC06dP17uOt7c3jh07huzsbHo4+gMqKyvRsWNH6OvrIz8/v3Xs/yfvKEhaP1ku2VBRPvl63yx3/PjxrKSkRNzuxIkTDABbuHBhg9dZunQpA8BOnjzZXF1XWAAa/fukCGgPD9JiMCrKp3BmzZqFpKQklJSUoKysDGlpaXB1dcXp06fh6ekp/l59zg7vpG2gIESaRXV1NXbv3o3w8HBs2rQJPB6vXhsqyqd4goOD4ezsjC5dukBTUxODBg3CiRMn4OTkhNTUVJw6dUreXSQtnALn9RFFUlhYiBkzZkgcGz9+PPbu3YsuXboA+LSifG9/TEX5WhYlJSXMnDkTycnJSElJgbu7+2ft8E7aBpoJEZmT5ZJNaWkpLfG0QKI3Fq9fvwbweTu8k7aBZkJE5t69vyNasnF2dkZycjJOnToFd3d3OfWOyMKVK1cAoEk7vJO2gWZCRC5ESzYAkJKSAoCK8imaO3fuQCAQ1DuenJyMqKgoqKqq4u9//zuAv0qd+/n5oby8HKtXr5Zov2bNGrx8+RJ+fn6tI+VYxhhjePTokby7ITU0EyJyI40lG3Nzc6SmpqKwsLDefSFa4pGtQ4cOITIyEqNHj4aJiQlUVVWRmZmJM2fOQElJCdHR0TAyMhK3DwwMxPHjxxEZGYmMjAzY2dmBx+MhLi4O1tbWCAwMlONoiLzQTIjIzceWbN72oaJ8AHDmzJl616eifLI1cuRI/O1vf8O9e/ewe/du/Pvf/wafz8fXX3+Ny5cvw8/PT6K9uro6kpKSwOVyce/ePaxbtw6ZmZngcrlISkqifePaKrk+pURaPT6fz16+fFnv+KVLl5iamhpTVVWVqJNCRfkIaVto2x4iUyEhIR9dsnn7HfO72/a8u2TT0LY9YWFhWLlypXjbnoqKChw4cACVlZWIj4/HyJEjm3vYhJBGoiBEZOrChQv4/fffcePGDRQVFaGqqgpdu3aFk5MTuFwuHBwc6n1OaWkpQkNDERsbK77X4+Pjg1WrVr03ySAmJgYbNmwAn8+HiooKHB0dsXr1atjb28t6iISQJqAgRAghRG4oMYEQQojcUBAi5DP88ccfGDt2LDp37owOHTqgZ8+emDJlCvLy8iTalZWVISAgAMbGxlBVVYWxsTECAgLEzzA1ZP/+/XBwcIC6ujp0dXXh5uaGa9euyXpIhMgFLccR8gkYY5g3bx62bNkCU1NTuLi4QFNTE/n5+bhw4QJiYmLg5OQEgGojEdIocsvLI0QB/frrrwwAW7BgAautra13/u10cKqNRMjH0UyIkEaqrKyEgYEBdHR0kJWVBWXl9284whiDgYEBysrKUFhYWG+vtO7du6Njx47Iy8sTb1Xz448/Ys2aNdi9ezd8fX0lrjd//nxER0cjPj6+wbIVhCgquidESCMlJCTgxYsX8Pb2Rl1dHY4cOYKIiAhER0dL1DgCqDYSIY1Fe8cR0kii5ABlZWUMGDAAWVlZ4nNKSkrgcrn417/+BYBqIxHSWDQTIqSRiouLAQDr1q2DlpYW0tPT8erVK1y8eBEWFhZYt24dNm3aBIBqIxHSWBSECGkkoVAIAFBRUcHRo0dhb28PDQ0NDBs2DLGxsVBSUsK6devk3EtCFAsFIUIaSTRLGThwILp37y5xzsrKCr169cKDBw8gEAioNhIhjURBiJBGsrS0BADo6Og0eF50vLKy8rNrI5WXl6OwsLBR7QlpDSgIEdJIot247969W+9cTU0N7t+/D3V1dejp6VFtJEIaiYIQIY1kamqKcePG4f79+9i2bZvEuYiICAgEAkyYMAHKysqfVc565syZUFZWRlhYmMSyHJ/Px549e2BqaopRo0bJdpCENDN6WJWQT/DgwQMMGTIExcXFcHd3R+/evZGRkYHExEQYGxsjLS1NnGJNtZEIaQS57tdAiALKzc1lM2bMYPr6+qx9+/bM0NCQLViwgBUVFdVrKxAIGJfLZYaGhuK2XC5XYlued+3bt48NHDiQdejQgWlra7Px48ez9PR0WQ6JELmhmRAhhBC5oXtChBBC5IaCECGEELmhIEQIIURuKAgRQgiRGwpCLUBUVBQ4HA5mzZrV4PmSkhJ07doVWlpayM3NbebeEUKI7FAQagH8/f0xbNgw7Ny5EydOnKh3fv78+SguLsb69ethZGQkhx4SQohsUIp2C/Hw4UP0798fmpqa4PP56NSpEwAgJiYG06ZNg6urK06dOiXnXhJCiHRREGpBfv/9dyxYsACTJ0/GgQMHUFBQACsrKzDGwOfz6+3cTAghio6CUAvCGMO4ceNw9uxZHDp0CLt27cKpU6ewd+9eTJs2Td7dI4QQqaMg1MLk5eWhb9++qK6uRnV1NSZMmIAjR47Iu1uEECITlJjQwhgaGuK7775DdXU1VFVVER0dLe8uEUKIzFAQamGKioqwdetWAEB1dTVOnjwp5x4RQojsUBBqYb799ls8f/4ckZGR6NSpE7hcLp4+fSrvbhFCiEz8D1q2PDufRi0DAAAAAElFTkSuQmCC", + "text/plain": [ + "<Figure size 640x480 with 1 Axes>" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "6335" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "vertex_data = load_data_ply(\"OGAgamodon.ply\")\n", + "my_scatter(vertex_data)\n", + "len(vertex_data)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "1583" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "mean_points = mean_points_func(vertex_data)\n", + "len(mean_points)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "vertex_data = load_data_npy(\"Agamodon_anguliceps_skeleton.npy\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Sorting points\n", + "* Skip this step if the data came from Blender\n", + "* When points are extracted using the skeletonize method they have to be sorted, else the interpolation and smoothing function wont work.\n", + "* This function takes a bit to load (a few seconds)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "e7acc894674c4805b8e14b3db4fec402", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "interactive(children=(FloatSlider(value=15.0, description='Threshold:', min=1.0, step=1.0), IntSlider(value=30…" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sorted_points = interactive_sorting_points(vertex_data)\n", + "sorted_points" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Interpolation and subsamling of points\n", + "* The function uses the simple linear interpolation\n", + "* Note that if the data came from Blender \"sorted_points.result\" needs to changed to \"mean_points\"" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "515c4d4883cb4bf78c8d423fa62692f6", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "interactive(children=(IntSlider(value=5, description='Interpolated Points:', max=20), IntSlider(value=20, desc…" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "interpolated_points = interactive_interp_points(sorted_points.result)\n", + "interpolated_points" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Smoothing\n", + "* Usually it works just fine with poly = 1, which means it fits a linear polynomie to the points.\n", + "* In generel the window length has to a relative higher value than the poly order, so keep that in mind.\n", + "* Also mind how much smoothing there is done (the window length) since the suture may loose its characteristics. " + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "91924aa3e6b143b79875fb68da2504f1", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "interactive(children=(IntSlider(value=10, description='Window_length:', max=50, min=1), IntSlider(value=1, des…" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "smoothed_points = interactive_smooth_points(interpolated_points.result)\n", + "smoothed_points" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Calculate radii\n", + "* This is always done on the smoothed curve " + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Radii shape: (329,)\n", + "equidistant_points shape: (329, 3)\n", + "min of radii: nan\n", + "max of radii: nan\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "c:\\Users\\aske_\\OneDrive - Danmarks Tekniske Universitet\\Skrivebord\\DTU\\4Semester\\Fagprojekt\\arbejde\\myHelperFunctions.py:31: RuntimeWarning: invalid value encountered in divide\n", + " toCircumsphereCenter = ((np.cross(p12_X_p13, p12) * np.dot(p13, p13)) + (np.cross(p13, p12_X_p13) * np.dot(p12, p12))) / (2 * np.dot(p12_X_p13, p12_X_p13))\n" + ] + } + ], + "source": [ + "centers, radii = calculate_radii(smoothed_points.result)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Finding the turning points\n", + "* This is just to get an idea of where the points are located this is NOT the actual amount of turning points" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "d3974a807b52461c880697f4980ac56b", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "interactive(children=(IntSlider(value=55, description='T:'), IntSlider(value=30, description='Elevation:', max…" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "thresholded_points = interactive_thresholed_points(smoothed_points.result, radii)\n", + "thresholded_points" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Find the actual amount of turning points" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "719311f87fbb476da46cdbca2aeb33b5", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "interactive(children=(IntSlider(value=60, description='T:'), IntSlider(value=13, description='Skip Indices:'),…" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "turning_points = interactive_turning_points(smoothed_points.result, radii)\n", + "turning_points" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " The distance between the start and end point is: 414.3270365155942 Voxels\n", + " The total length of the line is: 1140.332397731631 Voxels\n", + " The tortuosity of the line is: 2.7522519585533054\n", + " The volume of the bounding box is: 14687337.705951998 Voxels\n" + ] + }, + { + "data": { + "text/plain": [ + "(414.3270365155942, 1140.332397731631, 2.7522519585533054, 14687337.705951998)" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "curve_characteristics(smoothed_points.result)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Fagprojekt", + "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.10.14" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} -- GitLab