From 1e3cdbd3d9afd33f5002cdd01ac4881a9e3c1371 Mon Sep 17 00:00:00 2001 From: Tom Birdsong Date: Fri, 21 May 2021 11:41:48 -0400 Subject: [PATCH] ENH: Use hasi pipeline functions in example notebooks --- examples/MeshToMeshRegistration.ipynb | 548 ++++++---- examples/TemplateGenerationIterative.ipynb | 1106 +++++++++++++++----- src/hasi/hasi/align.py | 12 +- 3 files changed, 1153 insertions(+), 513 deletions(-) diff --git a/examples/MeshToMeshRegistration.ipynb b/examples/MeshToMeshRegistration.ipynb index e6f756a..c37ef32 100644 --- a/examples/MeshToMeshRegistration.ipynb +++ b/examples/MeshToMeshRegistration.ipynb @@ -6,7 +6,7 @@ "source": [ "# Mesh to Mesh Registration Example\n", "\n", - "ITK natively supports image-to-image registration, which is a common operation for medical images with symmetry. Another common method of storing 3D volumetric data is to represent volume surfaces with meshes. In this example we seek to register two meshes using various ITK metrics and optimization techniques.\n", + "ITK natively supports image-to-image registration, which is a common operation for medical images with symmetry. Another common method of storing 3D volumetric data is to represent volume surfaces as meshes. One way to get surface features from a sample mesh is to register a common atlas of correspondences and find related points on the surface of the sample mesh. In this example we seek to register two meshes using various ITK metrics and optimization techniques.\n", "\n", "Registration classes are defined in the Python `hasi` submodule and built on top of the ITK Python wrapping. The `MeanSquaresRegistrar` and `DiffeoRegistrar` classes apply registration techniques to images derived from mesh inputs, while the `PointSetEntropyRegistrar` aims to register meshes via point set entropy metrics. Mesh registration is carried out with each class in this notebook on sample bone femur mesh data downloaded to the `examples/Data` folder." ] @@ -24,13 +24,22 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import sys\n", + "!{sys.executable} -m pip install itk itkwidgets" + ] + }, + { + "cell_type": "code", + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Update sys.path to reference src/ modules\n", "import os\n", - "import sys\n", "import copy\n", "import importlib\n", "from urllib.request import urlretrieve\n", @@ -43,12 +52,16 @@ "module_path = os.path.abspath(os.path.join('..'))\n", "\n", "if module_path not in sys.path:\n", - " sys.path.append(module_path)" + " sys.path.append(module_path)\n", + " \n", + "# Ignore wrapping warnings from itkwidgets\n", + "import warnings\n", + "warnings.filterwarnings('ignore')" ] }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ @@ -58,7 +71,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ @@ -74,7 +87,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ @@ -89,7 +102,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 6, "metadata": { "scrolled": true }, @@ -103,40 +116,44 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Compare images with ITKWidgets\n", + "## Compare geometries with ITKWidgets\n", "\n", - "We can use `view`, `compare`, and `checkerboard` to inspect mesh and image data." + "We can use `view`, `compare`, and `checkerboard` to inspect mesh and image data. Comparing the two meshes, we see that they are generally similar but do not precisely align by default. Attempting to set sample correspondences from the template mesh with this default alignment would yield a poor description of the sample surface. Registration will align the surfaces so that the two bones better coincide." ] }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 7, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "3f97be0e263e437d8e479ea4c06bc1c9", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Viewer(geometries=[{'vtkClass': 'vtkPolyData', 'points': {'vtkClass': 'vtkPoints', 'numberOfComponents': 3, 'd…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ - "#view(geometries=[template_mesh,target_mesh])" + "view(geometries=[template_mesh,target_mesh])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Mesh-To-Image Conversion in the Base Class\n", - "Python registration classes inheriting from `MeshToMeshRegistrar` implement unique registration algorithms. The base class includes common definitions for 3D mesh-to-mesh registration, abstract methods, and a mesh-to-image conversion method.\n", + "## Mesh-To-Image Conversion\n", "\n", - "The `mesh_to_image` function takes in a 3D mesh and converts it into an ITK 3D image object. The spacing, origin, and size of the 3D image may be calculated from the minimum bounding box of the mesh or set from a reference image.\n", + "The `hasi` package provides Python functions for various stages of the HASI shape analysis pipeline in its `align` module. These include static methods for creating surface meshes from sample images, deforming an atlas to a sample, iteratively refining an atlas from a population, and more.\n", "\n", - "If a list of meshes is passed to `mesh_to_image` then the output images will each be set with the same size, spacing, and origin so that every mesh is fully contained within the image region." - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [], - "source": [ - "from src.hasi.hasi.meshtomeshregistrar import MeshToMeshRegistrar\n", - "registrar = MeshToMeshRegistrar()" + "The `align.mesh_to_image` function takes in a set of 3D meshes and outputs a set of ITK 3D images describing the meshes in a common space. The spacing, origin, and size of the images may be derived from the minimum common bounding box of the mesh or set from a reference image. The final space of the meshes is generated with a small empty buffer around the sample representation (5% in each direction)." ] }, { @@ -145,25 +162,16 @@ "metadata": {}, "outputs": [], "source": [ - "target_image, template_image = registrar.mesh_to_image([target_mesh, template_mesh])" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Comparing the meshes generated from the given meshes, we see that the two bone images are generally similar but do not exactly line up together. Registration will translate one mesh so that the two bones better coincide." + "from src.hasi.hasi.align import mesh_to_image" ] }, { "cell_type": "code", "execution_count": 9, - "metadata": { - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ - "#compare(template_image,target_image)" + "target_image, template_image = mesh_to_image([target_mesh, template_mesh])" ] }, { @@ -172,41 +180,50 @@ "metadata": { "scrolled": true }, - "outputs": [], - "source": [ - "#checkerboard(template_image, target_image, pattern=PATTERN_COUNT)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "86c55950b93f460ea8de289690969607", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "AppLayout(children=(HBox(children=(Label(value='Link:'), Checkbox(value=False, description='cmap'), Checkbox(v…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ - "## Point Set Downsizing in the Base Class\n", - "\n", - "The base `MeshToMeshRegistrar` class also provides decimation functionality for uniform random point set sampling from a given mesh, primarily used to improve performance for point set based registration of dense meshes. Note that uniform random sampling may not be suitable for all shapes, in which case application-specific resampling should be carried out externally prior to registration.\n", - "\n", - "For subsequent registration in this notebook we will rely on a template atlas with approximately 4,000 points which does not require decimation." + "compare(template_image,target_image)" ] }, { "cell_type": "code", "execution_count": 11, - "metadata": {}, - "outputs": [], - "source": [ - "target_points_reduced = \\\n", - " registrar.randomly_sample_mesh_points(mesh=target_mesh, sampling_rate=0.01)" - ] - }, - { - "cell_type": "code", - "execution_count": 12, "metadata": { - "scrolled": false + "scrolled": true }, - "outputs": [], + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "95b94e1d51ad4b7ebe414a9a8d135fe7", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "VBox(children=(Viewer(annotations=False, interpolation=False, rendered_image= list:\n", - " downsamples = list()\n", - " for image in image_list:\n", - " output_spacing = [spacing * ratio for spacing in itk.spacing(image)]\n", - " output_size = [int(size / ratio) for size in itk.size(image)]\n", - "\n", - " downsample = itk.resample_image_filter(Input=image,\n", - " Size=output_size,\n", - " OutputOrigin=itk.origin(image),\n", - " OutputSpacing=output_spacing,)\n", - " downsamples.append(downsample)\n", - " return downsamples" - ] - }, - { - "cell_type": "code", - "execution_count": 15, + "execution_count": 10, "metadata": {}, "outputs": [ { @@ -267,8 +189,8 @@ } ], "source": [ - "sparse_downsampled_images = downsample_images(images,SPARSE_DOWNSAMPLE_RATIO)\n", - "dense_downsampled_images = downsample_images(images,DENSE_DOWNSAMPLE_RATIO)\n", + "sparse_downsampled_images = align.downsample_images(images,SPARSE_DOWNSAMPLE_RATIO)\n", + "dense_downsampled_images = align.downsample_images(images,DENSE_DOWNSAMPLE_RATIO)\n", "\n", "print(itk.size(dense_downsampled_images[0]))\n", "print(itk.size(sparse_downsampled_images[0]))" @@ -276,7 +198,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 11, "metadata": { "scrolled": false }, @@ -284,7 +206,7 @@ { "data": { "application/vnd.jupyter.widget-view+json": { - "model_id": "3023646685ae428b92a13e8378fd0f8b", + "model_id": "6f25616a386b4871b8e0a328105932ea", "version_major": 2, "version_minor": 0 }, @@ -312,7 +234,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 12, "metadata": { "scrolled": true }, @@ -320,43 +242,25 @@ "source": [ "# Here each pixel interior to the femur has value \"1\" and exterior has value \"0\".\n", "# This may change for a different segmentation image.\n", - "FEMUR_OBJECT_PIXEL_VALUE = 1\n", - "\n", - "Dimension = itk.template(sparse_downsampled_images[0])[1][1]\n", - "MeshType = itk.Mesh[itk.F,Dimension]" + "FEMUR_OBJECT_PIXEL_VALUE = 1" ] }, { "cell_type": "code", - "execution_count": 16, - "metadata": {}, - "outputs": [], - "source": [ - "def generate_meshes(images:list, mesh_type=itk.Mesh[itk.F,3]) -> list:\n", - " meshes = list()\n", - " for image in images:\n", - " mesh = itk.binary_mask3_d_mesh_source(image, \n", - " object_value=1, \n", - " ttype=[type(images[0]),mesh_type])\n", - " meshes.append(mesh)\n", - " return meshes " - ] - }, - { - "cell_type": "code", - "execution_count": 17, + "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Average dense mesh points: 215290.5\n" + "Average dense mesh points: 215237.7857142857\n" ] } ], "source": [ - "dense_meshes = generate_meshes(dense_downsampled_images)\n", + "dense_meshes = align.binary_image_list_to_meshes(dense_downsampled_images,\n", + " object_pixel_value=FEMUR_OBJECT_PIXEL_VALUE)\n", "\n", "# Expect ~200K vertices\n", "print('Average dense mesh points: ' +\n", @@ -365,28 +269,29 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Average sparse mesh points: 4373.357142857143\n" + "Average sparse mesh points: 4371.5\n" ] } ], "source": [ - "sparse_meshes = generate_meshes(sparse_downsampled_images)\n", + "sparse_meshes = align.binary_image_list_to_meshes(sparse_downsampled_images,\n", + " object_pixel_value=FEMUR_OBJECT_PIXEL_VALUE)\n", "\n", - "# Expect ~5K vertices\n", + "# Expect ~200K vertices\n", "print('Average sparse mesh points: ' +\n", " str(sum([mesh.GetNumberOfPoints() for mesh in sparse_meshes]) / len(sparse_meshes)))" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 15, "metadata": { "scrolled": true }, @@ -402,59 +307,709 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 16, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Template itk::ImageToPointSetFilter,itk::PointSet>\n", + " already defined as \n", + " is redefined as {cl}\n", + "Warning: template already defined 'itk::ImageToPointSetFilter,itk::PointSet>'\n", + "Template itk::ImageToPointSetFilter,itk::PointSet>\n", + " already defined as \n", + " is redefined as {cl}\n", + "Warning: template already defined 'itk::ImageToPointSetFilter,itk::PointSet>'\n", + "Template itk::ImageToPointSetFilter,itk::PointSet>\n", + " already defined as \n", + " is redefined as {cl}\n", + "Warning: template already defined 'itk::ImageToPointSetFilter,itk::PointSet>'\n", + "Template itk::ImageToPointSetFilter,itk::PointSet>\n", + " already defined as \n", + " is redefined as {cl}\n", + "Warning: template already defined 'itk::ImageToPointSetFilter,itk::PointSet>'\n", + "Template itk::ImageToPointSetFilter,itk::PointSet>\n", + " already defined as \n", + " is redefined as {cl}\n", + "Warning: template already defined 'itk::ImageToPointSetFilter,itk::PointSet>'\n", + "Template itk::ImageToPointSetFilter,itk::PointSet>\n", + " already defined as \n", + " is redefined as {cl}\n", + "Warning: template already defined 'itk::ImageToPointSetFilter,itk::PointSet>'\n", + "Template itk::ImageToPointSetFilter,itk::PointSet>\n", + " already defined as \n", + " is redefined as {cl}\n", + "Warning: template already defined 'itk::ImageToPointSetFilter,itk::PointSet>'\n", + "Template itk::ImageToPointSetFilter,itk::PointSet>\n", + " already defined as \n", + " is redefined as {cl}\n", + "Warning: template already defined 'itk::ImageToPointSetFilter,itk::PointSet>'\n", + "Template itk::ImageToPointSetFilter,itk::PointSet>\n", + " already defined as \n", + " is redefined as {cl}\n", + "Warning: template already defined 'itk::ImageToPointSetFilter,itk::PointSet>'\n", + "Template itk::ImageToPointSetFilter,itk::PointSet>\n", + " already defined as \n", + " is redefined as {cl}\n", + "Warning: template already defined 'itk::ImageToPointSetFilter,itk::PointSet>'\n", + "Template itk::ImageToPointSetFilter,itk::PointSet>\n", + " already defined as \n", + " is redefined as {cl}\n", + "Warning: template already defined 'itk::ImageToPointSetFilter,itk::PointSet>'\n", + "Template itk::ImageToPointSetFilter,itk::PointSet>\n", + " already defined as \n", + " is redefined as {cl}\n", + "Warning: template already defined 'itk::ImageToPointSetFilter,itk::PointSet>'\n", + "Template itk::ImageToPointSetFilter,itk::PointSet>\n", + " already defined as \n", + " is redefined as {cl}\n", + "Warning: template already defined 'itk::ImageToPointSetFilter,itk::PointSet>'\n", + "Template itk::ImageToPointSetFilter,itk::PointSet>\n", + " already defined as \n", + " is redefined as {cl}\n", + "Warning: template already defined 'itk::ImageToPointSetFilter,itk::PointSet>'\n", + "Template itk::ImageToPointSetFilter,itk::PointSet>\n", + " already defined as \n", + " is redefined as {cl}\n", + "Warning: template already defined 'itk::ImageToPointSetFilter,itk::PointSet>'\n" + ] + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "48bf8486ae32422aa2080915a6b44da6", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Viewer(geometries=[{'vtkClass': 'vtkPolyData', 'points': {'vtkClass': 'vtkPoints', 'numberOfComponents': 3, 'd…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ "# visualize with itkwidgets\n", - "#view(geometries=sparse_meshes)" + "view(geometries=[mesh for mesh in sparse_meshes])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Define Iterative Registration\n", - "It is desired to find correspondence points between the template and each sample mesh. In order to get correspondence, two steps are employed:\n", + "### Register Sample for Correspondence\n", + "Our goal is to find correspondence points between the template and each sample mesh. In order to get correspondence, two steps are employed:\n", "- First, a copy of the template mesh is registered to the target sample;\n", "- Second, each mesh point is repositioned at its nearest neighbor on the target sample.\n", "\n", - "The result of this process is a full collection of meshes having the same number of points and correspondence between each point such that it represents the same approximate feature on each femur." + "The result of this process is a full collection of meshes having the same number of points and correspondence between each point such that it represents the same approximate feature on each femur. The cells below show an example of the template being registered to a single sample mesh. This process is repeated for every sample in the full pipeline." ] }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 17, "metadata": { "scrolled": true }, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Iteration: 0 Metric: 0.0\n", + "Iteration: 0 Metric: 0.0\n", + "Iteration: 0 Metric: 0.0\n", + "Iteration: 0 Metric: 0.0\n", + "Iteration: 0 Metric: 0.0\n", + "Iteration: 0 Metric: 0.7779061854942507\n", + "Iteration: 1 Metric: 0.20232000292614374\n", + "Iteration: 2 Metric: 0.13674457929331374\n", + "Iteration: 3 Metric: 0.12514542842076454\n", + "Iteration: 4 Metric: 0.122550291731051\n", + "Iteration: 5 Metric: 0.12061420415595286\n", + "Iteration: 6 Metric: 0.11841144247855505\n", + "Iteration: 7 Metric: 0.11589176266252948\n", + "Iteration: 8 Metric: 0.11325031330727194\n", + "Iteration: 9 Metric: 0.11063534863433273\n", + "Iteration: 10 Metric: 0.10807184441782834\n", + "Iteration: 11 Metric: 0.10561519395469256\n", + "Iteration: 12 Metric: 0.10323751692873938\n", + "Iteration: 13 Metric: 0.10090638903838738\n", + "Iteration: 14 Metric: 0.09868776219709642\n", + "Iteration: 15 Metric: 0.09658259228964905\n", + "Iteration: 16 Metric: 0.09458964512549776\n", + "Iteration: 17 Metric: 0.09270177175503178\n", + "Iteration: 18 Metric: 0.09092719465045711\n", + "Iteration: 19 Metric: 0.08921088449330158\n", + "Iteration: 20 Metric: 0.08751331057883362\n", + "Iteration: 21 Metric: 0.08595625062114805\n", + "Iteration: 22 Metric: 0.08448978369964545\n", + "Iteration: 23 Metric: 0.08306276567722867\n", + "Iteration: 24 Metric: 0.08171543662323255\n", + "Iteration: 25 Metric: 0.08043608000343186\n", + "Iteration: 26 Metric: 0.07922975383845902\n", + "Iteration: 27 Metric: 0.07806911113091398\n", + "Iteration: 28 Metric: 0.0769828791766969\n", + "Iteration: 29 Metric: 0.07592302414933781\n", + "Iteration: 30 Metric: 0.07494589189324528\n", + "Iteration: 31 Metric: 0.07396736060192158\n", + "Iteration: 32 Metric: 0.07305007692282282\n", + "Iteration: 33 Metric: 0.07217091810544289\n", + "Iteration: 34 Metric: 0.07137870580845175\n", + "Iteration: 35 Metric: 0.07062566206862178\n", + "Iteration: 36 Metric: 0.06990448486572416\n", + "Iteration: 37 Metric: 0.06922109639851765\n", + "Iteration: 38 Metric: 0.06853380709305265\n", + "Iteration: 39 Metric: 0.06788873436392066\n", + "Iteration: 40 Metric: 0.0672722282894907\n", + "Iteration: 41 Metric: 0.06669437953985918\n", + "Iteration: 42 Metric: 0.06615220704647713\n", + "Iteration: 43 Metric: 0.06561881914708118\n", + "Iteration: 44 Metric: 0.0651103655497439\n", + "Iteration: 45 Metric: 0.06463178159354313\n", + "Iteration: 46 Metric: 0.06419378162288007\n", + "Iteration: 47 Metric: 0.06378539931051558\n", + "Iteration: 48 Metric: 0.06339390676537018\n", + "Iteration: 49 Metric: 0.06299462478066936\n", + "Iteration: 50 Metric: 0.06259765910714457\n", + "Iteration: 51 Metric: 0.062223501410776164\n", + "Iteration: 52 Metric: 0.06186476948546618\n", + "Iteration: 53 Metric: 0.06151061722063784\n", + "Iteration: 54 Metric: 0.06118354489978483\n", + "Iteration: 55 Metric: 0.06087213131634895\n", + "Iteration: 56 Metric: 0.06057013719210121\n", + "Iteration: 57 Metric: 0.06027498100856217\n", + "Iteration: 58 Metric: 0.05999541178989893\n", + "Iteration: 59 Metric: 0.05972059722194064\n", + "Iteration: 60 Metric: 0.05944584050731294\n", + "Iteration: 61 Metric: 0.05917541864616805\n", + "Iteration: 62 Metric: 0.058908971747794144\n", + "Iteration: 63 Metric: 0.058655260566518784\n", + "Iteration: 64 Metric: 0.058415430581464645\n", + "Iteration: 65 Metric: 0.058178836264712354\n", + "Iteration: 66 Metric: 0.057937409843928145\n", + "Iteration: 67 Metric: 0.057697151558410986\n", + "Iteration: 68 Metric: 0.05747225797292385\n", + "Iteration: 69 Metric: 0.05727423976713176\n", + "Iteration: 70 Metric: 0.05709037866354649\n", + "Iteration: 71 Metric: 0.05691542628966276\n", + "Iteration: 72 Metric: 0.056746657780375326\n", + "Iteration: 73 Metric: 0.05657913987842849\n", + "Iteration: 74 Metric: 0.05641514413854398\n", + "Iteration: 75 Metric: 0.056251362421798484\n", + "Iteration: 76 Metric: 0.056090432428837184\n", + "Iteration: 77 Metric: 0.05593873693335362\n", + "Iteration: 78 Metric: 0.05579209223722084\n", + "Iteration: 79 Metric: 0.05565609855107295\n", + "Iteration: 80 Metric: 0.0555228740624395\n", + "Iteration: 81 Metric: 0.05538491139523176\n", + "Iteration: 82 Metric: 0.055250594503693316\n", + "Iteration: 83 Metric: 0.055121970213393746\n", + "Iteration: 84 Metric: 0.05500291647599293\n", + "Iteration: 85 Metric: 0.054891744576532966\n", + "Iteration: 86 Metric: 0.05478383310930921\n", + "Iteration: 87 Metric: 0.054681873426192164\n", + "Iteration: 88 Metric: 0.0545842727694528\n", + "Iteration: 89 Metric: 0.054489492014450656\n", + "Iteration: 90 Metric: 0.05438684625230607\n", + "Iteration: 91 Metric: 0.05429039237273384\n", + "Iteration: 92 Metric: 0.05420227931773402\n", + "Iteration: 93 Metric: 0.0541178113578612\n", + "Iteration: 94 Metric: 0.05402946310069399\n", + "Iteration: 95 Metric: 0.0539281294698308\n", + "Iteration: 96 Metric: 0.05383914528785153\n", + "Iteration: 97 Metric: 0.053758757355710056\n", + "Iteration: 98 Metric: 0.053674653015075825\n", + "Iteration: 99 Metric: 0.053584798945789866\n", + "Iteration: 100 Metric: 0.05349382055175694\n", + "Iteration: 101 Metric: 0.05341009320456643\n", + "Iteration: 102 Metric: 0.05332881383329767\n", + "Iteration: 103 Metric: 0.053256532548124506\n", + "Iteration: 104 Metric: 0.05319328769701666\n", + "Iteration: 105 Metric: 0.053131653243475706\n", + "Iteration: 106 Metric: 0.05307704756507718\n", + "Iteration: 107 Metric: 0.05302060581934514\n", + "Iteration: 108 Metric: 0.05295340984159919\n", + "Iteration: 109 Metric: 0.05289001870842689\n", + "Iteration: 110 Metric: 0.05282385749602994\n", + "Iteration: 111 Metric: 0.052765058546700064\n", + "Iteration: 112 Metric: 0.052697587776129824\n", + "Iteration: 113 Metric: 0.05263494878354366\n", + "Iteration: 114 Metric: 0.05256945628213697\n", + "Iteration: 115 Metric: 0.05250504946239966\n", + "Iteration: 116 Metric: 0.052450526747697526\n", + "Iteration: 117 Metric: 0.05239895469950103\n", + "Iteration: 118 Metric: 0.05235319450752982\n", + "Iteration: 119 Metric: 0.052312412068669024\n", + "Iteration: 120 Metric: 0.052268564087777705\n", + "Iteration: 121 Metric: 0.0522214106949559\n", + "Iteration: 122 Metric: 0.052173354704088175\n", + "Iteration: 123 Metric: 0.05212340589885686\n", + "Iteration: 124 Metric: 0.052078422799893984\n", + "Iteration: 125 Metric: 0.052031849859089076\n", + "Iteration: 126 Metric: 0.051988323846011\n", + "Iteration: 127 Metric: 0.051949387451224995\n", + "Iteration: 128 Metric: 0.05191534420597778\n", + "Iteration: 129 Metric: 0.05188183736711797\n", + "Iteration: 130 Metric: 0.051841031615754214\n", + "Iteration: 131 Metric: 0.05180304581125652\n", + "Iteration: 132 Metric: 0.05177053415449397\n", + "Iteration: 133 Metric: 0.05174146628941411\n", + "Iteration: 134 Metric: 0.05170791558239239\n", + "Iteration: 135 Metric: 0.05168041336473835\n", + "Iteration: 136 Metric: 0.051653407663689116\n", + "Iteration: 137 Metric: 0.05162798857126075\n", + "Iteration: 138 Metric: 0.051605632399284146\n", + "Iteration: 139 Metric: 0.05158414454641372\n", + "Iteration: 140 Metric: 0.05156124141667879\n", + "Iteration: 141 Metric: 0.051535505328572255\n", + "Iteration: 142 Metric: 0.05151744439212391\n", + "Iteration: 143 Metric: 0.0514928965718451\n", + "Iteration: 144 Metric: 0.05146978819894209\n", + "Iteration: 145 Metric: 0.051452394271285125\n", + "Iteration: 146 Metric: 0.05143333849263162\n", + "Iteration: 147 Metric: 0.05141683983858978\n", + "Iteration: 148 Metric: 0.05140006000618391\n", + "Iteration: 149 Metric: 0.051379632528070346\n", + "Iteration: 150 Metric: 0.051364163225486044\n", + "Iteration: 151 Metric: 0.051344148804672154\n", + "Iteration: 152 Metric: 0.0513209198248402\n", + "Iteration: 153 Metric: 0.051300204451203146\n", + "Iteration: 154 Metric: 0.051279847038341424\n", + "Iteration: 155 Metric: 0.05125897474166499\n", + "Iteration: 156 Metric: 0.05124039364195694\n", + "Iteration: 157 Metric: 0.051225711659785685\n", + "Iteration: 158 Metric: 0.051213232380085696\n", + "Iteration: 159 Metric: 0.05119908540041881\n", + "Iteration: 160 Metric: 0.05118816073223949\n", + "Iteration: 161 Metric: 0.05117790965210165\n", + "Iteration: 162 Metric: 0.051164989486057126\n", + "Iteration: 163 Metric: 0.05115162286003827\n", + "Iteration: 164 Metric: 0.051141654999141875\n", + "Iteration: 165 Metric: 0.05112915649279891\n", + "Iteration: 166 Metric: 0.05112134513222777\n", + "Iteration: 167 Metric: 0.051113475657624995\n", + "Iteration: 168 Metric: 0.05110545616598484\n", + "Iteration: 169 Metric: 0.05109767407846486\n", + "Iteration: 170 Metric: 0.051086077346439454\n", + "Iteration: 171 Metric: 0.05107590830530297\n", + "Iteration: 172 Metric: 0.05106466261954617\n", + "Iteration: 173 Metric: 0.05105306999550187\n", + "Iteration: 174 Metric: 0.05104524411816563\n", + "Iteration: 175 Metric: 0.051037250720081036\n", + "Iteration: 176 Metric: 0.051027554599183546\n", + "Iteration: 177 Metric: 0.051019245272718464\n", + "Iteration: 178 Metric: 0.0510103015813352\n", + "Iteration: 179 Metric: 0.05099999631723761\n", + "Iteration: 180 Metric: 0.05099747702422614\n", + "Iteration: 181 Metric: 0.05099025910556314\n", + "Iteration: 182 Metric: 0.05098246677270456\n", + "Iteration: 183 Metric: 0.05097646728829948\n", + "Iteration: 184 Metric: 0.0509680084125315\n", + "Iteration: 185 Metric: 0.050957756849324895\n", + "Iteration: 186 Metric: 0.05094888751493839\n", + "Iteration: 187 Metric: 0.05093697704947861\n", + "Iteration: 188 Metric: 0.05092846670001227\n", + "Iteration: 189 Metric: 0.050923196874734906\n", + "Iteration: 190 Metric: 0.0509168646827519\n", + "Iteration: 191 Metric: 0.050905025166283834\n", + "Iteration: 192 Metric: 0.05088763831250859\n", + "Iteration: 193 Metric: 0.05088060842698278\n", + "Iteration: 194 Metric: 0.05087075561956609\n", + "Iteration: 195 Metric: 0.050859520430219664\n", + "Iteration: 196 Metric: 0.05085411683245611\n", + "Iteration: 197 Metric: 0.05084410138767771\n", + "Iteration: 198 Metric: 0.05083532634158362\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Iteration: 199 Metric: 0.050829382864870494\n", + "Iteration: 200 Metric: 0.05082170914386962\n", + "Iteration: 201 Metric: 0.05081423866472783\n", + "Iteration: 202 Metric: 0.05080568023614978\n", + "Iteration: 203 Metric: 0.05079902262172653\n", + "Iteration: 204 Metric: 0.05079059523282873\n", + "Iteration: 205 Metric: 0.05078300527810398\n", + "Iteration: 206 Metric: 0.05077665585338589\n", + "Iteration: 207 Metric: 0.05076873757269327\n", + "Iteration: 208 Metric: 0.05076176368430392\n", + "Iteration: 209 Metric: 0.05076191191460448\n", + "Iteration: 210 Metric: 0.050759255781437584\n", + "Iteration: 211 Metric: 0.05075297240683363\n", + "Iteration: 212 Metric: 0.05074280619842757\n", + "Iteration: 213 Metric: 0.0507353605396458\n", + "Iteration: 214 Metric: 0.0507293774424416\n", + "Iteration: 215 Metric: 0.05072597680923839\n", + "Iteration: 216 Metric: 0.050725369103616885\n", + "Iteration: 217 Metric: 0.0507235730125808\n", + "Iteration: 218 Metric: 0.05071929349275945\n", + "Iteration: 219 Metric: 0.05071276157849241\n", + "Iteration: 220 Metric: 0.050704833611047065\n", + "Iteration: 221 Metric: 0.05070007184266171\n", + "Iteration: 222 Metric: 0.05069326759814266\n", + "Iteration: 223 Metric: 0.050687986995510546\n", + "Iteration: 224 Metric: 0.05068121272230936\n", + "Iteration: 225 Metric: 0.05067712628603081\n", + "Iteration: 226 Metric: 0.050671607763591114\n", + "Iteration: 227 Metric: 0.05066573637544863\n", + "Iteration: 228 Metric: 0.050661630055786895\n", + "Iteration: 229 Metric: 0.05065762856878695\n", + "Iteration: 230 Metric: 0.050653153150137\n", + "Iteration: 231 Metric: 0.05064987299222196\n", + "Iteration: 232 Metric: 0.05064315389881709\n", + "Iteration: 233 Metric: 0.05063952433441711\n", + "Iteration: 234 Metric: 0.05063466540346055\n", + "Iteration: 235 Metric: 0.05062943284705007\n", + "Iteration: 236 Metric: 0.05062446513810828\n", + "Iteration: 237 Metric: 0.0506187349444279\n", + "Iteration: 238 Metric: 0.05061558707153199\n", + "Iteration: 239 Metric: 0.05061220664007433\n", + "Iteration: 240 Metric: 0.050608453663761495\n", + "Iteration: 241 Metric: 0.05060405388796036\n", + "Iteration: 242 Metric: 0.05059907738048261\n", + "Iteration: 243 Metric: 0.05059543121583007\n", + "Iteration: 244 Metric: 0.05059200253936908\n", + "Iteration: 245 Metric: 0.05059161578220578\n", + "Iteration: 246 Metric: 0.05058699150323093\n", + "Iteration: 247 Metric: 0.05058291307147369\n", + "Iteration: 248 Metric: 0.050577805145974174\n", + "Iteration: 249 Metric: 0.05057217165358974\n", + "Iteration: 250 Metric: 0.050571062925596905\n", + "Iteration: 251 Metric: 0.05056496573044829\n", + "Iteration: 252 Metric: 0.050562110970127035\n", + "Iteration: 253 Metric: 0.05055728536395525\n", + "Iteration: 254 Metric: 0.05055235366588629\n", + "Iteration: 255 Metric: 0.05054737396130461\n", + "Iteration: 256 Metric: 0.05054440043124623\n", + "Iteration: 257 Metric: 0.05054188303247858\n", + "Iteration: 258 Metric: 0.05053881837291953\n", + "Iteration: 259 Metric: 0.05053234625570165\n", + "Iteration: 260 Metric: 0.05052936539212418\n", + "Iteration: 261 Metric: 0.05052522606071205\n", + "Iteration: 262 Metric: 0.050521127631690255\n", + "Iteration: 263 Metric: 0.05052246227319725\n", + "Iteration: 264 Metric: 0.050517423605369864\n", + "Iteration: 265 Metric: 0.05051349205007255\n", + "Iteration: 266 Metric: 0.05051099265176425\n", + "Iteration: 267 Metric: 0.050508973692926215\n", + "Iteration: 268 Metric: 0.05050551462038012\n", + "Iteration: 269 Metric: 0.050502524037637345\n", + "Iteration: 270 Metric: 0.050499477928054426\n", + "Iteration: 271 Metric: 0.05049629849772875\n", + "Iteration: 272 Metric: 0.05049760118512804\n", + "Iteration: 273 Metric: 0.050494756776694634\n", + "Iteration: 274 Metric: 0.05049055434412369\n", + "Iteration: 275 Metric: 0.05048739284126354\n", + "Iteration: 276 Metric: 0.050482198188618564\n", + "Iteration: 277 Metric: 0.05047929580808147\n", + "Iteration: 278 Metric: 0.0504751878945688\n", + "Iteration: 279 Metric: 0.05047077602731021\n", + "Iteration: 280 Metric: 0.050465679778863945\n", + "Iteration: 281 Metric: 0.050462028373489266\n", + "Iteration: 282 Metric: 0.05045805846907594\n", + "Iteration: 283 Metric: 0.050454642707376894\n", + "Iteration: 284 Metric: 0.050450230167405796\n", + "Iteration: 285 Metric: 0.050446779213573344\n", + "Iteration: 286 Metric: 0.050443380421668425\n", + "Iteration: 287 Metric: 0.05044124380681699\n", + "Iteration: 288 Metric: 0.0504394860849241\n", + "Iteration: 289 Metric: 0.05043804286975267\n", + "Iteration: 290 Metric: 0.050437299696372666\n", + "Iteration: 291 Metric: 0.050432846525755165\n", + "Iteration: 292 Metric: 0.05043123811402428\n", + "Iteration: 293 Metric: 0.05042600614502405\n", + "Iteration: 294 Metric: 0.05042245290959345\n", + "Iteration: 295 Metric: 0.050417389038791026\n", + "Iteration: 296 Metric: 0.05041356076936511\n", + "Iteration: 297 Metric: 0.05040913396334622\n", + "Iteration: 298 Metric: 0.05040420966281016\n", + "Iteration: 299 Metric: 0.05040278224274329\n", + "Iteration: 300 Metric: 0.050401768630244766\n", + "Iteration: 301 Metric: 0.05039872658215754\n", + "Iteration: 302 Metric: 0.05039429292407768\n", + "Iteration: 303 Metric: 0.05039157619714431\n", + "Iteration: 304 Metric: 0.05038822246984111\n", + "Iteration: 305 Metric: 0.050383375177383075\n", + "Iteration: 306 Metric: 0.0503805094230236\n", + "Iteration: 307 Metric: 0.05037832060080366\n", + "Iteration: 308 Metric: 0.050374129488421786\n", + "Iteration: 309 Metric: 0.05037167465089316\n", + "Iteration: 310 Metric: 0.050368475933972694\n", + "Iteration: 311 Metric: 0.05036502987680338\n", + "Iteration: 312 Metric: 0.05036328948124166\n", + "Iteration: 313 Metric: 0.05036304455490719\n", + "Iteration: 314 Metric: 0.050360530034586085\n", + "Iteration: 315 Metric: 0.05035837642959368\n", + "Iteration: 316 Metric: 0.05035485754947832\n", + "Iteration: 317 Metric: 0.05035390608912031\n", + "Iteration: 318 Metric: 0.05035212109171042\n", + "Iteration: 319 Metric: 0.05035103871150679\n", + "Iteration: 320 Metric: 0.05035011328891645\n", + "Iteration: 321 Metric: 0.05034859021203525\n", + "Iteration: 322 Metric: 0.050347311419667286\n", + "Iteration: 323 Metric: 0.05034482451091769\n", + "Iteration: 324 Metric: 0.050344502999895985\n", + "Iteration: 325 Metric: 0.05034167345627448\n", + "Iteration: 326 Metric: 0.05033954454776951\n", + "Iteration: 327 Metric: 0.050338444807378854\n", + "Iteration: 328 Metric: 0.05033333892514485\n", + "Iteration: 329 Metric: 0.05033061695326114\n", + "Iteration: 330 Metric: 0.050328536986903295\n", + "Iteration: 331 Metric: 0.050325284955912646\n", + "Iteration: 332 Metric: 0.05032113795012483\n", + "Iteration: 333 Metric: 0.050318440363230674\n", + "Iteration: 334 Metric: 0.050315692835759554\n", + "Iteration: 335 Metric: 0.050310221774776884\n", + "Iteration: 336 Metric: 0.050307172757623336\n", + "Iteration: 337 Metric: 0.050304335918295265\n", + "Iteration: 338 Metric: 0.050302998852957415\n", + "Iteration: 339 Metric: 0.050299817591392314\n", + "Iteration: 340 Metric: 0.05029759969848791\n", + "Iteration: 341 Metric: 0.05029510993358949\n", + "Iteration: 342 Metric: 0.05029122136606211\n", + "Iteration: 343 Metric: 0.05028921697481013\n", + "Iteration: 344 Metric: 0.05028675933115594\n", + "Iteration: 345 Metric: 0.050283812442764295\n", + "Iteration: 346 Metric: 0.050281927743758954\n", + "Iteration: 347 Metric: 0.05027902561494371\n", + "Iteration: 348 Metric: 0.05027841284200891\n", + "Iteration: 349 Metric: 0.050273692579973356\n", + "Iteration: 350 Metric: 0.05027127734164066\n", + "Iteration: 351 Metric: 0.05026854623742124\n", + "Iteration: 352 Metric: 0.05026598477111682\n", + "Iteration: 353 Metric: 0.05026377775027612\n", + "Iteration: 354 Metric: 0.05026080899514568\n", + "Iteration: 355 Metric: 0.050260336182794914\n", + "Iteration: 356 Metric: 0.05025836559014859\n", + "Iteration: 357 Metric: 0.050256441618175454\n", + "Iteration: 358 Metric: 0.05025404058669818\n", + "Iteration: 359 Metric: 0.050251581225139806\n", + "Iteration: 360 Metric: 0.050250054709434146\n", + "Iteration: 361 Metric: 0.05024844423755407\n", + "Iteration: 362 Metric: 0.05024696256938067\n", + "Iteration: 363 Metric: 0.050243679326953006\n", + "Iteration: 364 Metric: 0.05024205871054532\n", + "Iteration: 365 Metric: 0.05024065722054119\n", + "Iteration: 366 Metric: 0.05023820875244211\n", + "Iteration: 367 Metric: 0.05023613769674013\n", + "Iteration: 368 Metric: 0.05023555529977827\n", + "Iteration: 369 Metric: 0.050234782085916294\n", + "Iteration: 370 Metric: 0.05023372316288263\n", + "Iteration: 371 Metric: 0.05023041006198993\n", + "Iteration: 372 Metric: 0.0502283742613524\n", + "Iteration: 373 Metric: 0.05022586075701607\n", + "Iteration: 374 Metric: 0.050223304015505024\n", + "Iteration: 375 Metric: 0.05022193669254255\n", + "Iteration: 376 Metric: 0.05021741490878071\n", + "Iteration: 377 Metric: 0.05021548423049242\n", + "Iteration: 378 Metric: 0.05021318832704701\n", + "Iteration: 379 Metric: 0.05021125341836687\n", + "Iteration: 380 Metric: 0.050210314308740976\n", + "Iteration: 381 Metric: 0.050207382476498016\n", + "Iteration: 382 Metric: 0.050205047603897116\n", + "Iteration: 383 Metric: 0.05020377646884177\n", + "Iteration: 384 Metric: 0.050202725803133\n", + "Iteration: 385 Metric: 0.05020159735330152\n", + "Iteration: 386 Metric: 0.05019977700818119\n", + "Iteration: 387 Metric: 0.05019659475499572\n", + "Iteration: 388 Metric: 0.05019382863798433\n", + "Iteration: 389 Metric: 0.0501914220347015\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Iteration: 390 Metric: 0.05018866379176612\n", + "Iteration: 391 Metric: 0.0501866580990529\n", + "Iteration: 392 Metric: 0.05018544461556436\n", + "Iteration: 393 Metric: 0.050183944161473114\n", + "Iteration: 394 Metric: 0.05018252535440755\n", + "Iteration: 395 Metric: 0.050181952039632355\n", + "Iteration: 396 Metric: 0.05018052191351881\n", + "Iteration: 397 Metric: 0.05017906303736656\n", + "Iteration: 398 Metric: 0.050178842342465696\n", + "Iteration: 399 Metric: 0.05017926878989274\n", + "Iteration: 400 Metric: 0.0501761496165269\n", + "Iteration: 401 Metric: 0.05017516704134721\n", + "Iteration: 402 Metric: 0.05017371936068065\n", + "Iteration: 403 Metric: 0.05017193665910303\n", + "Iteration: 404 Metric: 0.050170784743674424\n", + "Iteration: 405 Metric: 0.050170514887143504\n", + "Iteration: 406 Metric: 0.050171765973031336\n", + "Iteration: 407 Metric: 0.050170869118313775\n", + "Iteration: 408 Metric: 0.05016937934734629\n", + "Iteration: 409 Metric: 0.0501673844188967\n", + "Iteration: 410 Metric: 0.050166100196711874\n", + "Iteration: 411 Metric: 0.05016514666894262\n", + "Iteration: 412 Metric: 0.05016612976804814\n", + "Iteration: 413 Metric: 0.05016357882674013\n", + "Iteration: 414 Metric: 0.050162055361429175\n", + "Iteration: 415 Metric: 0.05016021897741215\n", + "Iteration: 416 Metric: 0.050158974869405\n", + "Iteration: 417 Metric: 0.050158220964876724\n", + "Iteration: 418 Metric: 0.050157306725792966\n", + "Iteration: 419 Metric: 0.050155201441845036\n", + "Iteration: 420 Metric: 0.050160625074469035\n", + "Iteration: 421 Metric: 0.05015604999918623\n", + "Iteration: 422 Metric: 0.05015480240021365\n", + "Iteration: 423 Metric: 0.05015077919491364\n", + "Iteration: 424 Metric: 0.05014984714187697\n", + "Iteration: 425 Metric: 0.050147897342842265\n", + "Iteration: 426 Metric: 0.05014750205717864\n", + "Iteration: 427 Metric: 0.05014676747573058\n", + "Iteration: 428 Metric: 0.05014439829130655\n", + "Iteration: 429 Metric: 0.05014171150148085\n", + "Iteration: 430 Metric: 0.05013783652494085\n", + "Iteration: 431 Metric: 0.050136598736406475\n", + "Iteration: 432 Metric: 0.050136189352769554\n", + "Iteration: 433 Metric: 0.05013564703574842\n", + "Iteration: 434 Metric: 0.050137191495986944\n", + "Iteration: 435 Metric: 0.050135733190428446\n", + "Iteration: 436 Metric: 0.05013407258124941\n", + "Iteration: 437 Metric: 0.050131244740288856\n", + "Iteration: 438 Metric: 0.05013077696487938\n", + "Iteration: 439 Metric: 0.05012929740093765\n", + "Iteration: 440 Metric: 0.0501272224464538\n", + "Iteration: 441 Metric: 0.05012622208074147\n", + "Iteration: 442 Metric: 0.05012479969805858\n", + "Iteration: 443 Metric: 0.050124138938231866\n", + "Iteration: 444 Metric: 0.0501232227145544\n", + "Iteration: 445 Metric: 0.050122276680540256\n", + "Iteration: 446 Metric: 0.05012011797450067\n", + "Iteration: 447 Metric: 0.05011642153928056\n", + "Iteration: 448 Metric: 0.05011469844334742\n", + "Iteration: 449 Metric: 0.05011274779865784\n", + "Iteration: 450 Metric: 0.050110454865674454\n", + "Iteration: 451 Metric: 0.05010929522282516\n", + "Iteration: 452 Metric: 0.05010668016619564\n", + "Iteration: 453 Metric: 0.05010551075639925\n", + "Iteration: 454 Metric: 0.050103470218839916\n", + "Iteration: 455 Metric: 0.05010128360402523\n", + "Iteration: 456 Metric: 0.05009958402536302\n", + "Iteration: 457 Metric: 0.05009875797171815\n", + "Iteration: 458 Metric: 0.050099827762251886\n", + "Iteration: 459 Metric: 0.050098532874809996\n", + "Iteration: 460 Metric: 0.0500959213737748\n", + "Iteration: 461 Metric: 0.050094885953764406\n", + "Iteration: 462 Metric: 0.05009225552261694\n", + "Iteration: 463 Metric: 0.05009018220725623\n", + "Iteration: 464 Metric: 0.050089567358585066\n", + "Iteration: 465 Metric: 0.05008716836496022\n", + "Iteration: 466 Metric: 0.05008466078994798\n", + "Iteration: 467 Metric: 0.0500803773630747\n", + "Iteration: 468 Metric: 0.050075036843470366\n", + "Iteration: 469 Metric: 0.050071654360699755\n", + "Iteration: 470 Metric: 0.050070076813275695\n", + "Iteration: 471 Metric: 0.05006682416396176\n", + "Iteration: 472 Metric: 0.05006398214713756\n", + "Iteration: 473 Metric: 0.050063357436816475\n", + "Iteration: 474 Metric: 0.050062277529736375\n", + "Iteration: 475 Metric: 0.05006005400160558\n", + "Iteration: 476 Metric: 0.050057934210079764\n", + "Iteration: 477 Metric: 0.050056958310114794\n", + "Iteration: 478 Metric: 0.05005436748371105\n", + "Iteration: 479 Metric: 0.05005406283415361\n", + "Iteration: 480 Metric: 0.05005277249189518\n", + "Iteration: 481 Metric: 0.050051048384666064\n", + "Iteration: 482 Metric: 0.05005044755432356\n", + "Iteration: 483 Metric: 0.05004907424439143\n", + "Iteration: 484 Metric: 0.05004852043744791\n", + "Iteration: 485 Metric: 0.05004650240393613\n", + "Iteration: 486 Metric: 0.05004436681719005\n", + "Iteration: 487 Metric: 0.05004304873360667\n", + "Iteration: 488 Metric: 0.050041428661499326\n", + "Iteration: 489 Metric: 0.05004045188778732\n", + "Iteration: 490 Metric: 0.050039639222020954\n", + "Iteration: 491 Metric: 0.05003783421892475\n", + "Iteration: 492 Metric: 0.05003649428342877\n", + "Iteration: 493 Metric: 0.050036245476085435\n", + "Iteration: 494 Metric: 0.05003664703200972\n", + "Iteration: 495 Metric: 0.05003430333762224\n", + "Iteration: 496 Metric: 0.05003277583055426\n", + "Iteration: 497 Metric: 0.05003165744994026\n", + "Iteration: 498 Metric: 0.05003009200447867\n", + "Iteration: 499 Metric: 0.050030346055022885\n", + "Iteration: 500 Metric: 0.05003009200447867\n", + "Number of iterations run: 500\n", + "Final metric value: 0.05003009200447867\n", + "Final transform position: [0.8601696919525987, -0.06117567675444691, 0.1455632756241091, 0.024242047636822365, 0.9355363281830277, -0.03948992331825653, 0.024086151345177262, -0.003011293128523798, 0.9202711414886503, -0.9714223510362355, -0.3219046794482, -0.17710610686041697]\n", + "Optimizer scales: [54.022498598098714, 21.996100536727774, 18.232899837112292, 54.022498598098714, 21.996100536727774, 18.232899837112292, 54.02249859809877, 21.996100536727774, 18.23289983711267, 1.0000000000000018, 1.0000000000000018, 1.0000000000000018]\n", + "Optimizer learning rate: 1.0\n" + ] + } + ], "source": [ - "def register_template_to_sample(template_mesh, \n", - " sample_mesh,\n", - " learning_rate=1.0,\n", - " max_iterations=500):\n", - " \n", - " registrar = PointSetEntropyRegistrar()\n", - " metric = itk.EuclideanDistancePointSetToPointSetMetricv4[itk.PointSet[itk.F,3]].New()\n", - " transform = itk.Euler3DTransform[itk.D].New()\n", - " \n", - " # Make a deep copy of the template point set to register to the target\n", - " template_copy = itk.Mesh[itk.F,3].New()\n", - " for i in range(0, template_mesh.GetNumberOfPoints()):\n", - " template_copy.SetPoint(i, template_mesh.GetPoint(i))\n", - " template_copy.SetCells(template_mesh.GetCells())\n", - " \n", - " # Run registration and transform points to target\n", - " (transform, deformed_mesh) = registrar.register(template_mesh=template_copy,\n", - " target_mesh=sample_mesh,\n", - " metric=metric,\n", - " transform=transform,\n", - " learning_rate=learning_rate,\n", - " max_iterations=max_iterations,\n", - " resample_from_target=True,\n", - " verbose=False)\n", - " return deformed_mesh" + "registered_template = align.register_template_to_sample(sparse_meshes[1],\n", + " dense_meshes[2],\n", + " learning_rate=1.0,\n", + " minimum_convergence_value=-0.1,\n", + " max_iterations=500,\n", + " verbose=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "4cc752eb19bb4f1daa26feb935e83d38", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Viewer(geometries=[{'vtkClass': 'vtkPolyData', 'points': {'vtkClass': 'vtkPoints', 'numberOfComponents': 3, 'd…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "view(geometries=[registered_template, dense_meshes[2]])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Resample Template for Correspondence\n", + "\n", + "With the template mesh registered to the sample we can find correspondence points on the surface of the sample bone by getting the nearest sample point to every point on the template mesh surface." + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "e46e2901171c4e22a9a8d8ce2babbfe1", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Viewer(geometries=[{'vtkClass': 'vtkPolyData', 'points': {'vtkClass': 'vtkPoints', 'numberOfComponents': 3, 'd…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "sample_correspondence = align.resample_template_from_target(registered_template, dense_meshes[2])\n", + "view(geometries=[sample_correspondence,dense_meshes[2]])" ] }, { @@ -462,41 +1017,39 @@ "metadata": {}, "source": [ "### Define Procrustes Alignment Parameters\n", - "Once template meshes have been aligned to represent each input mesh with correspondence we can run Procrustes alignment and get out a mean mesh as the new template." + "Once template meshes have been aligned to represent each input mesh with correspondence we can run Procrustes alignment and get out a mean mesh as the new template. Here we demonstrate aligning two sample meshes as an example, but the full population is used in the atlas generation pipeline." ] }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 20, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "17d0ebec748b49f5a9cd17480fcb60db", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Viewer(geometries=[{'vtkClass': 'vtkPolyData', 'points': {'vtkClass': 'vtkPoints', 'numberOfComponents': 3, 'd…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ - "def align_procrustes(deformed_meshes, \n", - " template_index, \n", - " verbose=False,\n", - " convergence=0.08):\n", - " \n", - " procrustes_filter = itk.MeshProcrustesAlignFilter[type(deformed_meshes[0]), type(deformed_meshes[0])].New()\n", - " \n", - " procrustes_filter.SetUseInitialAverageOff()\n", - " procrustes_filter.SetUseNormalizationOff()\n", - " procrustes_filter.SetUseScalingOff()\n", - " procrustes_filter.SetConvergence(convergence) # Minimum threshold to exit alignment\n", - " \n", - " # Set mesh correspondence inputs\n", - " procrustes_filter.SetNumberOfInputs(len(deformed_meshes))\n", - " for i in range(0, len(deformed_meshes)):\n", - " procrustes_filter.SetInput(i, deformed_meshes[i])\n", - " \n", - " # Run filter\n", - " procrustes_filter.Update()\n", - " \n", - " if(verbose):\n", - " print(f'Alignment converged at {procrustes_filter.GetMeanPointsDifference()}')\n", - " \n", - " mean_result = procrustes_filter.GetMean()\n", - " mean_result.SetCells(deformed_meshes[template_index].GetCells())\n", - " return mean_result" + "second_registered_template = align.register_template_to_sample(sparse_meshes[1],\n", + " dense_meshes[3],\n", + " max_iterations=100)\n", + "second_sample_correspondence = align.resample_template_from_target(second_registered_template, dense_meshes[3])\n", + "\n", + "mean_correspondence = align.get_mean_correspondence_mesh([sample_correspondence,second_sample_correspondence],\n", + " convergence_threshold=1.0,\n", + " verbose=False)\n", + "view(geometries=[mean_correspondence])" ] }, { @@ -510,20 +1063,36 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 21, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Hausdorff distance between initial and updated mesh: 2.018958098608886\n" + ] + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "2211d8392c8f4787b1aff5a8769439ce", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Viewer(geometries=[{'vtkClass': 'vtkPolyData', 'points': {'vtkClass': 'vtkPoints', 'numberOfComponents': 3, 'd…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ - "def calculate_hausdorff_distance(mesh1, mesh2):\n", - " assert(mesh1.GetNumberOfPoints() == mesh2.GetNumberOfPoints())\n", - " max_dist = 0.0\n", - " \n", - " for pt_idx in range(mesh1.GetNumberOfPoints()):\n", - " pt1 = mesh1.GetPoint(pt_idx)\n", - " pt2 = mesh2.GetPoint(pt_idx)\n", - " dist = sum((pt1[dim] - pt2[dim]) ** 2 for dim in range(0,3)) ** 0.5\n", - " max_dist = max(max_dist, dist)\n", - " return max_dist" + "dist = align.get_pairwise_hausdorff_distance(mean_correspondence,sparse_meshes[1])\n", + "print(f'Hausdorff distance between initial and updated mesh: {dist}')\n", + "\n", + "view(geometries=[sparse_meshes[1],mean_correspondence])" ] }, { @@ -532,22 +1101,24 @@ "source": [ "### Run Iterative Refinement\n", "\n", - "Run registration and alignment on every mesh, choosing whether to iteratively update meshes to take advantage of mean alignment in each iteration or to ignore previous changes so that each alignment procedure is independent of the order of meshes in the list." + "Here we run the previously described registration, deformation, and alignment steps on every sample mesh, iterating the template as we go." ] }, { "cell_type": "code", - "execution_count": 23, + "execution_count": 22, "metadata": {}, "outputs": [], "source": [ - "# Fix the number of iterative refinements to run for each atlas.\n", - "# One iteration includes registration to each dense mesh, followed by\n", - "# subsequent Procrustes alignment of all correspondence meshes.\n", + "# Fix the number of iterative atlas refinements.\n", + "# For the mouse femur data set the atlas tends to converge\n", + "# after 3-5 iterations.\n", "NUM_ITERATIONS = 1\n", "\n", - "# Select indices of atlas templates to refine.\n", - "TEMPLATES_TO_ALIGN = [0]\n", + "# Select index of atlas template to refine.\n", + "# Initial template should be as close to a \"mean\" shape\n", + "# as possible, with no distinct surface outliers.\n", + "TEMPLATE_IDX = 1\n", "\n", "# Prepare directory to write out atlas iterations.\n", "# ex. 'Output/mean/0/901-L-femur-label.obj'\n", @@ -557,88 +1128,65 @@ }, { "cell_type": "code", - "execution_count": 24, - "metadata": { - "scrolled": true - }, + "execution_count": 23, + "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Now at iteration 0\n", - "Now at template mesh 0\n", - "Resampling from mesh 0\n", - "Resampling from mesh 1\n", - "Resampling from mesh 2\n", - "Resampling from mesh 3\n", - "Resampling from mesh 4\n", - "Resampling from mesh 5\n", - "Resampling from mesh 6\n", - "Resampling from mesh 7\n", - "Resampling from mesh 8\n", - "Resampling from mesh 9\n", - "Resampling from mesh 10\n", - "Resampling from mesh 11\n", - "Resampling from mesh 12\n", - "Resampling from mesh 13\n", - "Running alignment\n", - "Alignment converged at 0.0034345665327492492\n", - "Writing mean to Output/mean//0/901-R-femur-label.obj\n", - "Elapsed: 137.288982629776\n", - "Mesh 0 Hausdorff distance from previous iteration: 2.228128978646963\n" + "Time elapsed for iteration 0: 164.92673254013062\n", + "Writing to Output/mean//0/902-R-femur-label.obj\n" ] } ], "source": [ - "aligned_templates = dict()\n", - "for iteration in range(0,NUM_ITERATIONS):\n", - " print(f'Now at iteration {iteration}')\n", + "template_mesh = sparse_meshes[TEMPLATE_IDX]\n", + "for iteration in range(0, NUM_ITERATIONS):\n", + " starttime = time.time()\n", + " updated_template = align.refine_template_from_population(template_mesh=template_mesh,\n", + " target_meshes=dense_meshes,\n", + " registration_iterations=200)\n", + " endtime = time.time()\n", + " print(f'Time elapsed for iteration {iteration}: {endtime - starttime}')\n", " \n", - " for template_mesh_index in TEMPLATES_TO_ALIGN:\n", - " starttime = time.time()\n", - " print(f'Now at template mesh {template_mesh_index}')\n", - " template_mesh = sparse_meshes[template_mesh_index]\n", - "\n", - " deformed_templates = list()\n", - " for sample_index in range(0,len(dense_meshes)):\n", - " print(f'Resampling from mesh {sample_index}')\n", - " deformed_template = register_template_to_sample(template_mesh, \n", - " dense_meshes[sample_index],\n", - " max_iterations=50)\n", - " deformed_templates.append(deformed_template)\n", - "\n", - " print('Running alignment')\n", - " mesh_result = align_procrustes(deformed_templates, template_mesh_index, verbose=True)\n", - "\n", - " # Save intermediate results to disk\n", - " # These meshes are very small so this is not a significant expense (~400 KB/mesh)\n", - " output_path = f'{MEAN_OUTPUT_FOLDER}/{iteration}/{MESH_FILENAMES[template_mesh_index]}'\n", - " print(f'Writing mean to {output_path}')\n", - " #itk.meshwrite(mesh_result, output_path)\n", - " \n", - " # Optionally update current mesh in place for use in subsequent alignments\n", - " aligned_templates[template_mesh_index] = mesh_result\n", - " \n", - " endtime = time.time()\n", - " print(f'Elapsed: {endtime - starttime}')\n", + " output_path = f'{MEAN_OUTPUT_FOLDER}/{iteration}/{MESH_FILENAMES[TEMPLATE_IDX]}'\n", + " print(f'Writing to {output_path}')\n", + " itk.meshwrite(updated_template, output_path)\n", " \n", - " # Optionally update templates only between iterations\n", - " for template_mesh_index in TEMPLATES_TO_ALIGN:\n", - " dist = calculate_hausdorff_distance(sparse_meshes[template_mesh_index],\n", - " aligned_templates[template_mesh_index])\n", - " print(f'Mesh {template_mesh_index} Hausdorff distance from previous iteration: {dist}')\n", - " sparse_meshes[template_mesh_index] = aligned_templates[template_mesh_index]" + " template_mesh = updated_template" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The final template may be further examined via point distance metrics to determine fitness as an atlas. Here we will simply view the result." ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 24, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "6145223f116848f794ae1da95a839b5b", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Viewer(geometries=[{'vtkClass': 'vtkPolyData', 'points': {'vtkClass': 'vtkPoints', 'numberOfComponents': 3, 'd…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ - "# visualize template\n", - "# view(geometries=sparse_meshes)" + "# Visualize \n", + "view(geometries=[template_mesh])" ] } ], diff --git a/src/hasi/hasi/align.py b/src/hasi/hasi/align.py index c027441..192c4d3 100644 --- a/src/hasi/hasi/align.py +++ b/src/hasi/hasi/align.py @@ -200,14 +200,16 @@ def binary_image_list_to_meshes(images:list, mesh_pixel_type:type=itk.F, object_ # Register template mesh to sample mesh def register_template_to_sample(template_mesh:itk.Mesh, sample_mesh:itk.Mesh, - learning_rate=1.0, - max_iterations=2000, - verbose=False): + transform=None, + learning_rate:float=1.0, + max_iterations:int=2000, + minimum_convergence_value:float=-1.0, + convergence_window_size:int=1, + verbose=False) -> itk.Mesh: from .pointsetentropyregistrar import PointSetEntropyRegistrar registrar = PointSetEntropyRegistrar(verbose=verbose) metric = itk.EuclideanDistancePointSetToPointSetMetricv4[itk.PointSet[itk.F,3]].New() - transform = itk.Euler3DTransform[itk.D].New() # Make a deep copy of the template point set to resample from the target template_copy = deep_copy_mesh_vertices(template_mesh) @@ -218,6 +220,8 @@ def register_template_to_sample(template_mesh:itk.Mesh, metric=metric, transform=transform, learning_rate=learning_rate, + minimum_convergence_value=minimum_convergence_value, + convergence_window_size=convergence_window_size, max_iterations=max_iterations) return deformed_mesh