diff --git a/examples/TemplateGenerationIterative.ipynb b/examples/TemplateGenerationIterative.ipynb
index d0e5634..6fd3530 100644
--- a/examples/TemplateGenerationIterative.ipynb
+++ b/examples/TemplateGenerationIterative.ipynb
@@ -6,7 +6,9 @@
    "source": [
     "# Iterative Template Generation\n",
     "\n",
-    "This notebook exemplifies one way in which a template mesh atlas can be generated from a collection of segmented binary images. Each binary image of a mouse femur is downsampled to reduce pixel density prior to applying marching cubes to generate a mesh from the binary image. One arbitrary mesh is selected as the template and then registered to and resampled from each original mesh to get a full set of meshes with correspondence points. The meshes are then groupwise registered via procrustes alignment and the mean mesh is taken as the new template. This process is repeated for a fixed number of iterations to get a template mesh atlas that represents the average case of all meshes."
+    "This notebook exemplifies one way in which a template mesh atlas can be generated from a collection of segmented binary images. Each binary image of a mouse femur is downsampled to reduce pixel density prior to applying marching cubes to generate a mesh from the binary image. One arbitrary mesh is selected as the template and then registered to and resampled from each original mesh to get a full set of meshes with correspondence points. The meshes are then groupwise registered via procrustes alignment and the mean mesh is taken as the new template. This process is repeated for a fixed number of iterations to get a template mesh atlas that represents the average case of all meshes.\n",
+    "\n",
+    "This pipeline uses the [ITKShape](https://github.com/slicersalt/ITKShape) module for shape analysis."
    ]
   },
   {
@@ -17,6 +19,7 @@
    "source": [
     "import os\n",
     "import glob\n",
+    "import time\n",
     "\n",
     "import numpy as np\n",
     "import itk\n",
@@ -37,7 +40,7 @@
    "source": [
     "### Read images\n",
     "\n",
-    "Input images represent the results of automatic binary segmentations of mouse femur data. Each image contains only the femur object and represents a different possible region in space."
+    "Input images represent the results of automatic binary segmentations of mouse femur data. Each image contains only the femur object and potentially represents a different region in space."
    ]
   },
   {
@@ -46,38 +49,41 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# TODO download femurs if not locally available\n",
     "IMAGE_FOLDER = 'Input/femurs/'\n",
-    "MESH_OUTPUT_FOLDER = 'Output/femurs/'\n",
-    "MEAN_OUTPUT_FOLDER = 'Output/mean/'"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 3,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "for folder in [IMAGE_FOLDER, MESH_OUTPUT_FOLDER, MEAN_OUTPUT_FOLDER]:\n",
+    "DENSE_MESH_OUTPUT_FOLDER = 'Output/femurs/'\n",
+    "TEMPLATE_OUTPUT_FOLDER = 'Output/templates/'\n",
+    "MEAN_OUTPUT_FOLDER = 'Output/mean/'\n",
+    "\n",
+    "for folder in [IMAGE_FOLDER, DENSE_MESH_OUTPUT_FOLDER, TEMPLATE_OUTPUT_FOLDER, MEAN_OUTPUT_FOLDER]:\n",
     "    os.makedirs(folder, exist_ok=True)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 3,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "['Input/femurs\\\\901-R-femur-label.nrrd', 'Input/femurs\\\\902-R-femur-label.nrrd', 'Input/femurs\\\\906-R-femur-label.nrrd', 'Input/femurs\\\\907-R-femur-label.nrrd', 'Input/femurs\\\\908-R-femur-label.nrrd', 'Input/femurs\\\\915-R-femur-label.nrrd', 'Input/femurs\\\\916-R-femur-label.nrrd', 'Input/femurs\\\\917-R-femur-label.nrrd', 'Input/femurs\\\\918-R-femur-label.nrrd', 'Input/femurs\\\\F9-3wk-01-R-femur-label.nrrd', 'Input/femurs\\\\F9-3wk-02-R-femur-label.nrrd', 'Input/femurs\\\\F9-3wk-03-R-femur-label.nrrd', 'Input/femurs\\\\F9-8wk-01-R-femur-label.nrrd', 'Input/femurs\\\\F9-8wk-02-R-femur-label.nrrd']\n"
+     ]
+    }
+   ],
    "source": [
-    "input_paths = glob.glob(IMAGE_FOLDER + '*')\n",
-    "assert(len(input_paths) == 28)\n",
+    "# Get healthy femur segmentation binary images at \n",
+    "# https://data.kitware.com/#collection/5dcc6691e3566bda4b802172/folder/5e0b8d6baf2e2eed35c326f7\n",
+    "\n",
+    "input_paths = glob.glob(IMAGE_FOLDER + '*-R-*')\n",
+    "assert(len(input_paths) == 14)\n",
     "\n",
-    "# FIXME remove 901-R for now because it is misaligned\n",
-    "del input_paths[1]"
+    "print(input_paths)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 4,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -86,12 +92,11 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 5,
    "metadata": {},
    "outputs": [],
    "source": [
     "images = list()\n",
-    "\n",
     "for path in input_paths:\n",
     "    images.append(itk.imread(path, itk.UC))"
    ]
@@ -101,12 +106,13 @@
    "metadata": {},
    "source": [
     "### Paste images into same space\n",
-    "For viewing convenience."
+    "\n",
+    "Here we standardize physical space across the femur images. This is primarily intended to assist in viewing convenience with `itkwidgets` which expects a standard viewing region, but could also be helpful to standardize output from subsequent image downsampling and mesh conversion operations."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 6,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -119,14 +125,14 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": 7,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "itkSize3 ([1392, 983, 1247])\n"
+      "itkSize3 ([1279, 954, 1039])\n"
      ]
     }
    ],
@@ -141,7 +147,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": 8,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -157,19 +163,19 @@
     "    new_image.SetSpacing(spacing)\n",
     "    new_image.Allocate()\n",
     "    \n",
-    "    filter = itk.PasteImageFilter[type(orig_image)].New()\n",
-    "    filter.SetSourceImage(orig_image)\n",
-    "    filter.SetSourceRegion(orig_image.GetLargestPossibleRegion())\n",
-    "    filter.SetDestinationImage(new_image)\n",
-    "    filter.SetDestinationIndex([0,0,0])\n",
-    "    filter.Update()\n",
+    "    paste_filter = itk.PasteImageFilter[type(orig_image)].New()\n",
+    "    paste_filter.SetSourceImage(orig_image)\n",
+    "    paste_filter.SetSourceRegion(orig_image.GetLargestPossibleRegion())\n",
+    "    paste_filter.SetDestinationImage(new_image)\n",
+    "    paste_filter.SetDestinationIndex([0,0,0])\n",
+    "    paste_filter.Update()\n",
     "    \n",
-    "    return filter.GetOutput()"
+    "    return paste_filter.GetOutput()"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 10,
+   "execution_count": 9,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -180,13 +186,15 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 11,
-   "metadata": {},
+   "execution_count": 10,
+   "metadata": {
+    "scrolled": true
+   },
    "outputs": [
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "ee821e422e2746df91da3d8c03bcb11f",
+       "model_id": "c24adb830a56485292265239f102b44e",
        "version_major": 2,
        "version_minor": 0
       },
@@ -199,6 +207,7 @@
     }
    ],
    "source": [
+    "# View a 3D image with itkwidgets\n",
     "view(images[1])"
    ]
   },
@@ -207,47 +216,67 @@
    "metadata": {},
    "source": [
     "### Downsample images\n",
-    "The marching cubes algorithm returns a mesh with vertex density related to the pixel density of the original image. In this case marching cubes on the default images would produce meshes of approximately 800,000 points each, but the template mesh atlas is desired to contain only approximately 5,000 points. Each image is downsampled in order to yield a less dense mesh output."
+    "\n",
+    "The marching cubes algorithm returns a mesh with vertex density related to the pixel density of the original image. Here we downsample each image twice, once to get a \"dense\" image retaining most information density and a second time to get a \"sparse\" image more easily applied to get correspondence points.\n",
+    "\n",
+    "Meshes generated later from the dense images have approximately 600,000 vertices each while meshes generated from the sparse images have approximately 4,000 vertices. We will use the dense meshes to sample feature information and iteratively refine a the atlas to generalize the shape population. We can select a single sparse mesh to act as the initial atlas or carry out iterative refinement on multiple sparse meshes and compare to determine which result \"best\" reflects the population."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 12,
+   "execution_count": 11,
    "metadata": {},
    "outputs": [],
    "source": [
-    "DOWNSAMPLE_RATIO = 14"
+    "SPARSE_DOWNSAMPLE_RATIO = 14\n",
+    "DENSE_DOWNSAMPLE_RATIO = 2"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": 12,
    "metadata": {},
    "outputs": [],
    "source": [
-    "downsampled_images = list()"
+    "def downsample_images(image_list, ratio) -> 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": 14,
+   "execution_count": 15,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "itkSize3 ([639, 477, 519])\n",
+      "itkSize3 ([91, 68, 74])\n"
+     ]
+    }
+   ],
    "source": [
-    "for image in images:\n",
-    "    output_spacing = [spacing * DOWNSAMPLE_RATIO for spacing in itk.spacing(image)]\n",
-    "    output_size = [int(size / DOWNSAMPLE_RATIO) for size in itk.size(image)]\n",
+    "sparse_downsampled_images = downsample_images(images,SPARSE_DOWNSAMPLE_RATIO)\n",
+    "dense_downsampled_images = downsample_images(images,DENSE_DOWNSAMPLE_RATIO)\n",
     "\n",
-    "    downsample = itk.resample_image_filter(Input=image,\n",
-    "                                           Size=output_size,\n",
-    "                                           OutputOrigin=itk.origin(image),\n",
-    "                                           OutputSpacing=output_spacing,)\n",
-    "    downsampled_images.append(downsample)"
+    "print(itk.size(dense_downsampled_images[0]))\n",
+    "print(itk.size(sparse_downsampled_images[0]))"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 15,
+   "execution_count": 16,
    "metadata": {
     "scrolled": false
    },
@@ -255,7 +284,7 @@
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "271e62f8d21146f2acc8818f912bb61e",
+       "model_id": "3023646685ae428b92a13e8378fd0f8b",
        "version_major": 2,
        "version_minor": 0
       },
@@ -268,7 +297,7 @@
     }
    ],
    "source": [
-    "view(downsample)"
+    "view(dense_downsampled_images[0])"
    ]
   },
   {
@@ -276,22 +305,41 @@
    "metadata": {},
    "source": [
     "### Generate Meshes\n",
-    "The `itk.BinaryMask3DMeshSource` class makes use of the Marching Cubes algorithm to generate a mesh from a given object. Each binary image here uses the value '1' to indicate the femur is present at a pixel and '0' to indicate the femur is not present. Marching Cubes rapidly fills the femur space and generates surfaces at pixel region boundaries."
+    "The `itk.BinaryMask3DMeshSource` class makes use of the Marching Cubes algorithm to generate a mesh from a given object. Each binary image here uses the value '1' to indicate the femur is present at a pixel and '0' to indicate the femur is not present. Marching Cubes rapidly fills the femur space and generates surfaces at pixel region boundaries.\n",
+    "\n",
+    "Note that it may be useful to visually examine intermediate results. In the case where a mesh is not well aligned with others it is useful to correct the transformation in an external mesh editor."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": 17,
    "metadata": {
     "scrolled": true
    },
    "outputs": [],
    "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(downsampled_images[0])[1][1]\n",
-    "MeshType = itk.Mesh[itk.F,Dimension]\n",
-    "MeshSourceType = itk.BinaryMask3DMeshSource[type(downsampled_images[0]), MeshType]"
+    "Dimension = itk.template(sparse_downsampled_images[0])[1][1]\n",
+    "MeshType = itk.Mesh[itk.F,Dimension]"
+   ]
+  },
+  {
+   "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    "
    ]
   },
   {
@@ -303,73 +351,63 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Average mesh points: 4822.518518518518\n"
+      "Average dense mesh points: 215290.5\n"
      ]
     }
    ],
    "source": [
-    "meshes = list()\n",
+    "dense_meshes = generate_meshes(dense_downsampled_images)\n",
     "\n",
-    "for image in downsampled_images:\n",
-    "    mesh_source = MeshSourceType.New()\n",
-    "    mesh_source.SetObjectValue(FEMUR_OBJECT_PIXEL_VALUE)\n",
-    "    mesh_source.SetInput(image)\n",
-    "    mesh_source.Update()\n",
-    "    \n",
-    "    mesh_output = mesh_source.GetOutput()\n",
-    "    meshes.append(mesh_output)\n",
-    "    \n",
-    "print('Average mesh points: ' +\n",
-    "      str(sum([mesh.GetNumberOfPoints() for mesh in meshes]) / len(meshes)))"
+    "# Expect ~200K vertices\n",
+    "print('Average dense mesh points: ' +\n",
+    "      str(sum([mesh.GetNumberOfPoints() for mesh in dense_meshes]) / len(dense_meshes)))"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": 18,
    "metadata": {},
-   "outputs": [],
-   "source": [
-    "# Write out each mesh to disk\n",
-    "for i in range(0,len(meshes)):\n",
-    "    itk.meshwrite(meshes[i], f'{MESH_OUTPUT_FOLDER}{MESH_FILENAMES[i]}')"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 19,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# TODO visualize with itkwidgets\n",
-    "# view(geometries=meshes)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Average sparse mesh points: 4373.357142857143\n"
+     ]
+    }
+   ],
    "source": [
-    "### Select Template Mesh\n",
-    "It is now necessary to select a mesh that will act as the 'standard' for correspondence point updates from this point forward. The first mesh in the list is arbitrarily selected as this standard template."
+    "sparse_meshes = generate_meshes(sparse_downsampled_images)\n",
+    "\n",
+    "# Expect ~5K 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": 20,
-   "metadata": {},
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
    "outputs": [],
    "source": [
-    "TEMPLATE_MESH_INDEX = 0\n",
-    "template_mesh = meshes[TEMPLATE_MESH_INDEX]"
+    "# Write out each mesh to disk\n",
+    "for i in range(0,len(dense_meshes)):\n",
+    "    itk.meshwrite(dense_meshes[i], f'{DENSE_MESH_OUTPUT_FOLDER}{MESH_FILENAMES[i]}')\n",
+    "\n",
+    "for i in range(0,len(sparse_meshes)):\n",
+    "    itk.meshwrite(sparse_meshes[i], f'{TEMPLATE_OUTPUT_FOLDER}{MESH_FILENAMES[i]}')"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 21,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
-    "# FIXME visualize template\n",
-    "# view(geometries=[template_mesh])"
+    "# visualize with itkwidgets\n",
+    "#view(geometries=sparse_meshes)"
    ]
   },
   {
@@ -379,31 +417,14 @@
     "### 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",
     "- First, a copy of the template mesh is registered to the target sample;\n",
-    "- Second, each mesh point is updated to its nearest neighbor on 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."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 22,
-   "metadata": {
-    "scrolled": true
-   },
-   "outputs": [],
-   "source": [
-    "#FIXME set mesh1 to better inital transform for registration\n",
-    "#transform = itk.Euler3DTransform[itk.D].New()\n",
-    "#params = transform.GetParameters()\n",
-    "#params.SetElement(0,3.14)\n",
-    "#params.SetElement(1,3.14)\n",
-    "#params.SetElement(2,3.14 / 2)\n",
-    "#transform.SetParameters(params)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 23,
+   "execution_count": 20,
    "metadata": {
     "scrolled": true
    },
@@ -412,19 +433,19 @@
     "def register_template_to_sample(template_mesh, \n",
     "                                sample_mesh,\n",
     "                                learning_rate=1.0,\n",
-    "                                max_iterations=2000):\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 resample from the target\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 resample from target\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",
@@ -441,12 +462,12 @@
    "metadata": {},
    "source": [
     "### Define Procrustes Alignment Parameters\n",
-    "Now that template meshes have been aligned to represent each input mesh with correspondence, 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."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 24,
+   "execution_count": 21,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -478,6 +499,33 @@
     "    return mean_result"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Compute Hausdorff Distance\n",
+    "\n",
+    "We can calculate the Hausdorff distance from the current to previous mesh atlas at each iterative refinement to quantify the amount of change between iterations. In this case the meshes are in correspondence so we get the largest Euclidean distance between any pair of correspondence points."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 22,
+   "metadata": {},
+   "outputs": [],
+   "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"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -489,15 +537,19 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 26,
+   "execution_count": 23,
    "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",
     "NUM_ITERATIONS = 1\n",
-    "UPDATE_CONTINUOUSLY = True\n",
-    "aligned_templates = list()\n",
     "\n",
-    "# Write out mean meshes by iteration\n",
+    "# Select indices of atlas templates to refine.\n",
+    "TEMPLATES_TO_ALIGN = [0]\n",
+    "\n",
+    "# Prepare directory to write out atlas iterations.\n",
     "# ex. 'Output/mean/0/901-L-femur-label.obj'\n",
     "for i in range(0,NUM_ITERATIONS):\n",
     "    os.makedirs(MEAN_OUTPUT_FOLDER + str(i) + '\\\\', exist_ok=True)"
@@ -505,7 +557,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 29,
+   "execution_count": 24,
    "metadata": {
     "scrolled": true
    },
@@ -516,402 +568,6 @@
      "text": [
       "Now at iteration 0\n",
       "Now at template 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",
-      "Resampling from mesh 14\n",
-      "Resampling from mesh 15\n",
-      "Resampling from mesh 16\n",
-      "Resampling from mesh 17\n",
-      "Resampling from mesh 18\n",
-      "Resampling from mesh 19\n",
-      "Resampling from mesh 20\n",
-      "Resampling from mesh 21\n",
-      "Resampling from mesh 22\n",
-      "Resampling from mesh 23\n",
-      "Resampling from mesh 24\n",
-      "Resampling from mesh 25\n",
-      "Resampling from mesh 26\n",
-      "Running alignment\n",
-      "Alignment converged at 0.0013539076418325581\n",
-      "Writing mean to Output/mean//0/901-L-femur-label.obj\n",
-      "Now at template mesh 1\n",
-      "Resampling from mesh 0\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",
-      "Resampling from mesh 14\n",
-      "Resampling from mesh 15\n",
-      "Resampling from mesh 16\n",
-      "Resampling from mesh 17\n",
-      "Resampling from mesh 18\n",
-      "Resampling from mesh 19\n",
-      "Resampling from mesh 20\n",
-      "Resampling from mesh 21\n",
-      "Resampling from mesh 22\n",
-      "Resampling from mesh 23\n",
-      "Resampling from mesh 24\n",
-      "Resampling from mesh 25\n",
-      "Resampling from mesh 26\n",
-      "Running alignment\n",
-      "Alignment converged at 0.0021126325521304205\n",
-      "Writing mean to Output/mean//0/902-L-femur-label.obj\n",
-      "Now at template mesh 2\n",
-      "Resampling from mesh 0\n",
-      "Resampling from mesh 1\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",
-      "Resampling from mesh 14\n",
-      "Resampling from mesh 15\n",
-      "Resampling from mesh 16\n",
-      "Resampling from mesh 17\n",
-      "Resampling from mesh 18\n",
-      "Resampling from mesh 19\n",
-      "Resampling from mesh 20\n",
-      "Resampling from mesh 21\n",
-      "Resampling from mesh 22\n",
-      "Resampling from mesh 23\n",
-      "Resampling from mesh 24\n",
-      "Resampling from mesh 25\n",
-      "Resampling from mesh 26\n",
-      "Running alignment\n",
-      "Alignment converged at 0.0013403246222149828\n",
-      "Writing mean to Output/mean//0/902-R-femur-label.obj\n",
-      "Now at template mesh 3\n",
-      "Resampling from mesh 0\n",
-      "Resampling from mesh 1\n",
-      "Resampling from mesh 2\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",
-      "Resampling from mesh 14\n",
-      "Resampling from mesh 15\n",
-      "Resampling from mesh 16\n",
-      "Resampling from mesh 17\n",
-      "Resampling from mesh 18\n",
-      "Resampling from mesh 19\n",
-      "Resampling from mesh 20\n",
-      "Resampling from mesh 21\n",
-      "Resampling from mesh 22\n",
-      "Resampling from mesh 23\n",
-      "Resampling from mesh 24\n",
-      "Resampling from mesh 25\n",
-      "Resampling from mesh 26\n",
-      "Running alignment\n",
-      "Alignment converged at 0.0010776738666368856\n",
-      "Writing mean to Output/mean//0/906-L-femur-label.obj\n",
-      "Now at template mesh 4\n",
-      "Resampling from mesh 0\n",
-      "Resampling from mesh 1\n",
-      "Resampling from mesh 2\n",
-      "Resampling from mesh 3\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",
-      "Resampling from mesh 14\n",
-      "Resampling from mesh 15\n",
-      "Resampling from mesh 16\n",
-      "Resampling from mesh 17\n",
-      "Resampling from mesh 18\n",
-      "Resampling from mesh 19\n",
-      "Resampling from mesh 20\n",
-      "Resampling from mesh 21\n",
-      "Resampling from mesh 22\n",
-      "Resampling from mesh 23\n",
-      "Resampling from mesh 24\n",
-      "Resampling from mesh 25\n",
-      "Resampling from mesh 26\n",
-      "Running alignment\n",
-      "Alignment converged at 0.0009368941564440909\n",
-      "Writing mean to Output/mean//0/906-R-femur-label.obj\n",
-      "Now at template mesh 5\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 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",
-      "Resampling from mesh 14\n",
-      "Resampling from mesh 15\n",
-      "Resampling from mesh 16\n",
-      "Resampling from mesh 17\n",
-      "Resampling from mesh 18\n",
-      "Resampling from mesh 19\n",
-      "Resampling from mesh 20\n",
-      "Resampling from mesh 21\n",
-      "Resampling from mesh 22\n",
-      "Resampling from mesh 23\n",
-      "Resampling from mesh 24\n",
-      "Resampling from mesh 25\n",
-      "Resampling from mesh 26\n",
-      "Running alignment\n",
-      "Alignment converged at 0.01771205656103327\n",
-      "Writing mean to Output/mean//0/907-L-femur-label.obj\n",
-      "Now at template mesh 6\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 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",
-      "Resampling from mesh 14\n",
-      "Resampling from mesh 15\n",
-      "Resampling from mesh 16\n",
-      "Resampling from mesh 17\n",
-      "Resampling from mesh 18\n",
-      "Resampling from mesh 19\n",
-      "Resampling from mesh 20\n",
-      "Resampling from mesh 21\n",
-      "Resampling from mesh 22\n",
-      "Resampling from mesh 23\n",
-      "Resampling from mesh 24\n",
-      "Resampling from mesh 25\n",
-      "Resampling from mesh 26\n",
-      "Running alignment\n",
-      "Alignment converged at 0.0015883326966217473\n",
-      "Writing mean to Output/mean//0/907-R-femur-label.obj\n",
-      "Now at template mesh 7\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 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",
-      "Resampling from mesh 14\n",
-      "Resampling from mesh 15\n",
-      "Resampling from mesh 16\n",
-      "Resampling from mesh 17\n",
-      "Resampling from mesh 18\n",
-      "Resampling from mesh 19\n",
-      "Resampling from mesh 20\n",
-      "Resampling from mesh 21\n",
-      "Resampling from mesh 22\n",
-      "Resampling from mesh 23\n",
-      "Resampling from mesh 24\n",
-      "Resampling from mesh 25\n",
-      "Resampling from mesh 26\n",
-      "Running alignment\n",
-      "Alignment converged at 0.0019276822157436202\n",
-      "Writing mean to Output/mean//0/908-L-femur-label.obj\n",
-      "Now at template mesh 8\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 9\n",
-      "Resampling from mesh 10\n",
-      "Resampling from mesh 11\n",
-      "Resampling from mesh 12\n",
-      "Resampling from mesh 13\n",
-      "Resampling from mesh 14\n",
-      "Resampling from mesh 15\n",
-      "Resampling from mesh 16\n",
-      "Resampling from mesh 17\n",
-      "Resampling from mesh 18\n",
-      "Resampling from mesh 19\n",
-      "Resampling from mesh 20\n",
-      "Resampling from mesh 21\n",
-      "Resampling from mesh 22\n",
-      "Resampling from mesh 23\n",
-      "Resampling from mesh 24\n",
-      "Resampling from mesh 25\n",
-      "Resampling from mesh 26\n",
-      "Running alignment\n",
-      "Alignment converged at 0.0008468561680195561\n",
-      "Writing mean to Output/mean//0/908-R-femur-label.obj\n",
-      "Now at template mesh 9\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 10\n",
-      "Resampling from mesh 11\n",
-      "Resampling from mesh 12\n",
-      "Resampling from mesh 13\n",
-      "Resampling from mesh 14\n",
-      "Resampling from mesh 15\n",
-      "Resampling from mesh 16\n",
-      "Resampling from mesh 17\n",
-      "Resampling from mesh 18\n",
-      "Resampling from mesh 19\n",
-      "Resampling from mesh 20\n",
-      "Resampling from mesh 21\n",
-      "Resampling from mesh 22\n",
-      "Resampling from mesh 23\n",
-      "Resampling from mesh 24\n",
-      "Resampling from mesh 25\n",
-      "Resampling from mesh 26\n",
-      "Running alignment\n",
-      "Alignment converged at 0.0014259736570786917\n",
-      "Writing mean to Output/mean//0/915-L-femur-label.obj\n",
-      "Now at template mesh 10\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 11\n",
-      "Resampling from mesh 12\n",
-      "Resampling from mesh 13\n",
-      "Resampling from mesh 14\n",
-      "Resampling from mesh 15\n",
-      "Resampling from mesh 16\n",
-      "Resampling from mesh 17\n",
-      "Resampling from mesh 18\n",
-      "Resampling from mesh 19\n",
-      "Resampling from mesh 20\n",
-      "Resampling from mesh 21\n",
-      "Resampling from mesh 22\n",
-      "Resampling from mesh 23\n",
-      "Resampling from mesh 24\n",
-      "Resampling from mesh 25\n",
-      "Resampling from mesh 26\n"
-     ]
-    },
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Running alignment\n",
-      "Alignment converged at 0.0007549640420944888\n",
-      "Writing mean to Output/mean//0/915-R-femur-label.obj\n",
-      "Now at template mesh 11\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 12\n",
-      "Resampling from mesh 13\n",
-      "Resampling from mesh 14\n",
-      "Resampling from mesh 15\n",
-      "Resampling from mesh 16\n",
-      "Resampling from mesh 17\n",
-      "Resampling from mesh 18\n",
-      "Resampling from mesh 19\n",
-      "Resampling from mesh 20\n",
-      "Resampling from mesh 21\n",
-      "Resampling from mesh 22\n",
-      "Resampling from mesh 23\n",
-      "Resampling from mesh 24\n",
-      "Resampling from mesh 25\n",
-      "Resampling from mesh 26\n",
-      "Running alignment\n",
-      "Alignment converged at 0.0008808021989129835\n",
-      "Writing mean to Output/mean//0/916-L-femur-label.obj\n",
-      "Now at template mesh 12\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 13\n",
-      "Resampling from mesh 14\n",
-      "Resampling from mesh 15\n",
-      "Resampling from mesh 16\n",
-      "Resampling from mesh 17\n",
-      "Resampling from mesh 18\n",
-      "Resampling from mesh 19\n",
-      "Resampling from mesh 20\n",
-      "Resampling from mesh 21\n",
-      "Resampling from mesh 22\n",
-      "Resampling from mesh 23\n",
-      "Resampling from mesh 24\n",
-      "Resampling from mesh 25\n",
-      "Resampling from mesh 26\n",
-      "Running alignment\n",
-      "Alignment converged at 0.0004644934773637309\n",
-      "Writing mean to Output/mean//0/916-R-femur-label.obj\n",
-      "Now at template mesh 13\n",
       "Resampling from mesh 0\n",
       "Resampling from mesh 1\n",
       "Resampling from mesh 2\n",
@@ -925,67 +581,32 @@
       "Resampling from mesh 10\n",
       "Resampling from mesh 11\n",
       "Resampling from mesh 12\n",
-      "Resampling from mesh 14\n",
-      "Resampling from mesh 15\n",
-      "Resampling from mesh 16\n",
-      "Resampling from mesh 17\n",
-      "Resampling from mesh 18\n",
-      "Resampling from mesh 19\n",
-      "Resampling from mesh 20\n",
-      "Resampling from mesh 21\n",
-      "Resampling from mesh 22\n",
-      "Resampling from mesh 23\n",
-      "Resampling from mesh 24\n",
-      "Resampling from mesh 25\n",
-      "Resampling from mesh 26\n",
+      "Resampling from mesh 13\n",
       "Running alignment\n",
-      "Alignment converged at 0.003948452450873146\n",
-      "Writing mean to Output/mean//0/917-L-femur-label.obj\n",
-      "Now at template mesh 14\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"
-     ]
-    },
-    {
-     "ename": "KeyboardInterrupt",
-     "evalue": "",
-     "output_type": "error",
-     "traceback": [
-      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
-      "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
-      "\u001b[1;32m<ipython-input-29-28689507c5ee>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     12\u001b[0m             \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     13\u001b[0m                 \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34mf'Resampling from mesh {sample_index}'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 14\u001b[1;33m                 \u001b[0mdeformed_template\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mregister_template_to_sample\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtemplate_mesh\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmeshes\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0msample_index\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     15\u001b[0m                 \u001b[0mdeformed_templates\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdeformed_template\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     16\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
-      "\u001b[1;32m<ipython-input-23-2a81ba112558>\u001b[0m in \u001b[0;36mregister_template_to_sample\u001b[1;34m(template_mesh, sample_mesh, learning_rate, max_iterations)\u001b[0m\n\u001b[0;32m     15\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     16\u001b[0m     \u001b[1;31m# Run registration and resample from target\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 17\u001b[1;33m     (transform, deformed_mesh) = registrar.register(template_mesh=template_copy,\n\u001b[0m\u001b[0;32m     18\u001b[0m                                                     \u001b[0mtarget_mesh\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0msample_mesh\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     19\u001b[0m                                                     \u001b[0mmetric\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mmetric\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
-      "\u001b[1;32m~\\source\\repos\\HASI\\src\\hasi\\hasi\\pointsetentropyregistrar.py\u001b[0m in \u001b[0;36mregister\u001b[1;34m(self, template_mesh, template_point_set, target_mesh, target_point_set, filepath, verbose, metric, transform, resample_rate, max_iterations, learning_rate, resample_from_target)\u001b[0m\n\u001b[0;32m    143\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    144\u001b[0m         \u001b[1;32mif\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mresample_from_target\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 145\u001b[1;33m             \u001b[0mregistered_template_mesh\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mresample_template_from_target\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mregistered_template_mesh\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtarget_mesh\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    146\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    147\u001b[0m         \u001b[1;31m# Write out\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
-      "\u001b[1;32m~\\source\\repos\\HASI\\src\\hasi\\hasi\\meshtomeshregistrar.py\u001b[0m in \u001b[0;36mresample_template_from_target\u001b[1;34m(template_mesh, target_mesh)\u001b[0m\n\u001b[0;32m    119\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    120\u001b[0m         \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtemplate_mesh\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mGetNumberOfPoints\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 121\u001b[1;33m             \u001b[0mpt\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mtemplate_mesh\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mGetPoint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    122\u001b[0m             \u001b[0mind\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mloc\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mFindClosestPoint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mpt\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    123\u001b[0m             \u001b[0mwarped_pts\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mSetPoint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtarget_pts\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mGetPoint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mind\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
-      "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
+      "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"
      ]
     }
    ],
    "source": [
+    "aligned_templates = dict()\n",
     "for iteration in range(0,NUM_ITERATIONS):\n",
     "    print(f'Now at iteration {iteration}')\n",
     "    \n",
-    "    for template_mesh_index in range(0,len(meshes)):\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 = meshes[template_mesh_index]\n",
+    "        template_mesh = sparse_meshes[template_mesh_index]\n",
     "\n",
     "        deformed_templates = list()\n",
-    "        for sample_index in range(0,len(meshes)):\n",
-    "            if sample_index == template_mesh_index:\n",
-    "                deformed_templates.append(template_mesh)\n",
-    "            else:\n",
-    "                print(f'Resampling from mesh {sample_index}')\n",
-    "                deformed_template = register_template_to_sample(template_mesh, meshes[sample_index])\n",
-    "                deformed_templates.append(deformed_template)\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",
@@ -994,17 +615,20 @@
     "        # 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",
+    "        #itk.meshwrite(mesh_result, output_path)\n",
     "        \n",
     "        # Optionally update current mesh in place for use in subsequent alignments\n",
-    "        if UPDATE_IN_PLACE:\n",
-    "            meshes[template_mesh_index] = mesh_result\n",
-    "        else:\n",
-    "            aligned_templates.append(mesh_result)\n",
+    "        aligned_templates[template_mesh_index] = mesh_result\n",
+    "        \n",
+    "        endtime = time.time()\n",
+    "        print(f'Elapsed: {endtime - starttime}')\n",
     "    \n",
     "    # Optionally update templates only between iterations\n",
-    "    if not UPDATE_CONTINUOUSLY:\n",
-    "        meshes = aligned_templates"
+    "    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]"
    ]
   },
   {
@@ -1013,16 +637,16 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# FIXME visualize template\n",
-    "# view(geometries=[meshes])"
+    "# visualize template\n",
+    "# view(geometries=sparse_meshes)"
    ]
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python (dev)",
+   "display_name": "Python 3",
    "language": "python",
-   "name": "venv-itkdev"
+   "name": "python3"
   },
   "language_info": {
    "codemirror_mode": {
@@ -1034,7 +658,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.8"
+   "version": "3.8.5"
   }
  },
  "nbformat": 4,
diff --git a/src/hasi/hasi/pointsetentropyregistrar.py b/src/hasi/hasi/pointsetentropyregistrar.py
index 00ec4cb..c4bcf59 100644
--- a/src/hasi/hasi/pointsetentropyregistrar.py
+++ b/src/hasi/hasi/pointsetentropyregistrar.py
@@ -20,7 +20,6 @@
 
 import logging
 import os
-import pytest
 
 import itk
 from .meshtomeshregistrar import MeshToMeshRegistrar
diff --git a/src/hasi/pyproject.toml b/src/hasi/pyproject.toml
index 373a217..444fec5 100644
--- a/src/hasi/pyproject.toml
+++ b/src/hasi/pyproject.toml
@@ -8,6 +8,7 @@ author = "Kitware Medical"
 author-email = "matt.mccormick@kitware.com"
 home-page = "https://github.com/KitwareMedical/HASI"
 requires = [
+    "itk-shape>=0.2.0",
     "itk-hasi"
 ]
 requires-python=">=3.7"