Repository: KitwareMedical/HASI (language: python)

Issue:
Add Python script for atlas generation
Summation of steps described in [atlas generation notebook](https://github.com/KitwareMedical/HASI/blob/master/examples/TemplateGenerationIterative.ipynb).


Discussion / hints:

Fix patch:
diff --git a/src/hasi/hasi/align.py b/src/hasi/hasi/align.py
new file mode 100644
index 0000000..a575295
--- /dev/null
+++ b/src/hasi/hasi/align.py
@@ -0,0 +1,306 @@
+#!/usr/bin/env python3
+
+# Copyright NumFOCUS
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#        http://www.apache.org/licenses/LICENSE-2.0.txt
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+# Purpose: Python functions for mesh registration and alignment
+#           in HASI pipeline
+
+import itk
+
+
+# Convert itk.Mesh object to 3-dimensional itk.Image object
+def mesh_to_image(meshes:list, reference_image:itk.Image=None) -> itk.Image :
+    # Allow single mesh as input
+    if type(meshes) != list:
+        meshes = [meshes]
+
+    mesh_type = type(meshes[0])
+    pixel_type, dimension = itk.template(mesh_type)[1]
+    image_type = itk.Image[pixel_type, dimension]
+
+    mesh_to_image_filter_type = itk.TriangleMeshToBinaryImageFilter[mesh_type,image_type]
+    images = list()
+
+    if not reference_image:
+        # Set bounds to largest region encompassing all meshes
+        # Bounds format: [x_min x_max y_min y_max z_min z_max]
+        bounds = meshes[0].GetBoundingBox().GetBounds()
+        for mesh in meshes[1:]:
+            mesh_bounds = mesh.GetBoundingBox().GetBounds()
+            for dim in range(0, dimension):
+                bounds[dim*2] = min(bounds[dim*2], mesh_bounds[dim*2])
+                bounds[(dim*2)+1] = max(bounds[(dim*2)+1], mesh_bounds[(dim*2)+1])
+
+        # Calculate spacing, origin, and size with 5 pixel buffer around mesh
+        spacing = ((bounds[1] - bounds[0]) / 90,
+                    (bounds[3] - bounds[2]) / 90,
+                    (bounds[5] - bounds[4]) / 90)
+        origin = (bounds[0] - 5 * spacing[0],
+                    bounds[2] - 5 * spacing[1],
+                    bounds[4] - 5 * spacing[2])
+        size = (100,100,100)
+        direction = itk.Matrix[itk.D,dimension,dimension].GetIdentity()
+    else:
+        # Use given parameters
+        origin = reference_image.GetOrigin()
+        spacing = reference_image.GetSpacing()
+        size = reference_image.GetLargestPossibleRegion().GetSize()
+        direction = reference_image.GetDirection()
+
+    # Generate image for each mesh
+    for mesh in meshes:
+        mesh_to_image_filter = mesh_to_image_filter_type.New(Input=mesh,
+                                                             Origin=origin,
+                                                             Spacing=spacing,
+                                                             Size=size,
+                                                             Direction=direction)
+        mesh_to_image_filter.Update()
+
+        distance = itk.signed_maurer_distance_map_image_filter(mesh_to_image_filter.GetOutput())
+        images.append(distance)
+
+    return images[0] if len(images) == 1 else images
+
+
+# Build KdTree for operations referencing the target mesh
+def build_kd_tree(target_mesh:itk.Mesh, bucket_size=16) -> itk.KdTree:
+    _, Dimension = itk.template(target_mesh)[1]
+
+    # Set samples from target mesh points
+    sample = itk.ListSample[itk.Vector[itk.F, Dimension]].New()
+
+    # Populate ListSample with mesh points
+    sample.Resize(target_mesh.GetNumberOfPoints())
+    for i in range(0, target_mesh.GetNumberOfPoints()):
+        sample.SetMeasurementVector(i, target_mesh.GetPoint(i))
+
+    # Build kd-tree
+    tree_generator = itk.KdTreeGenerator[type(sample)].New(Sample=sample, BucketSize=bucket_size)
+    tree_generator.Update()
+    tree = tree_generator.GetOutput()
+
+    # Iteratively update template points to nearest neighbor in target point set
+    return tree, sample
+
+
+# Deep copy itk.Mesh vertices to avoid updating template in place
+def deep_copy_mesh_vertices(mesh:itk.Mesh) -> itk.Mesh:
+    mesh_copy = type(mesh).New()
+
+    # Deep copy vertices
+    for i in range(mesh.GetNumberOfPoints()):
+        mesh_copy.SetPoint(i, mesh.GetPoint(i))
+
+    # Shallow copy cells
+    mesh_copy.SetCells(mesh.GetCells())
+
+    return mesh_copy
+
+
+# Align each template vertex to the closest target vertex
+def resample_template_from_target(template_mesh:itk.Mesh, target_mesh:itk.Mesh) -> itk.Mesh:
+    # Reference the KdTree to set each template vertex to its nearest neighbor in the target mesh
+    _, Dimension = itk.template(template_mesh)[1]
+    measurement = itk.Vector[itk.F,Dimension]()
+
+    # Use available wrapping for itk.VectorContainer
+    if (itk.ULL, itk.ULL) in itk.VectorContainer.keys():
+        neighbors = itk.VectorContainer[itk.ULL, itk.ULL].New()
+    else:
+        neighbors = itk.VectorContainer[itk.UL, itk.UL].New()
+
+    # Create a KdTree from the target
+    tree, _ = build_kd_tree(target_mesh=target_mesh, bucket_size=16)
+
+    template_copy = deep_copy_mesh_vertices(template_mesh)
+    for i in range(0, template_mesh.GetNumberOfPoints()):
+        tree.Search(template_mesh.GetPoint(i), 1, neighbors)
+        template_copy.SetPoint(i, target_mesh.GetPoint(neighbors.GetElement(0)))
+
+    return template_copy
+
+
+# Resize binary images of equal spacing to occupy the largest common region
+def paste_to_common_space(images:list) -> list:
+    image_type = type(images[0])
+    pixel_type, dimension = itk.template(images[0])[1]
+
+    resized_images = list()
+
+    # Verify spacing is equivalent
+    SPACING_TOLERANCE = 1e-7
+    assert(all([itk.spacing(images[idx])[dim] - itk.spacing(images[0])[dim] < SPACING_TOLERANCE
+                for dim in range(dimension)
+                for idx in range(1,len(images))]))
+
+    # Get largest common region
+    max_size = itk.size(images[0])
+    for image in images:
+        max_size = [max(max_size[i], itk.size(image)[i]) for i in range(len(max_size))]
+
+    # Define paste region
+    for image in images:
+        region = itk.ImageRegion[dimension]()
+        region.SetSize(max_size)
+        region.SetIndex([0] * dimension)
+        new_image = type(images[0]).New(regions=region, spacing=image.GetSpacing())
+        new_image.Allocate()
+
+        resized_image = itk.paste_image_filter(source_image=image,
+                                               source_region=image.GetLargestPossibleRegion(),
+                                               destination_image=new_image,
+                                               destination_index=[0] * dimension,
+                                               ttype=type(image))
+        resized_images.append(resized_image)
+    return resized_images
+
+
+# Downsample images by a given ratio
+def downsample_images(images:list, ratio:float) -> list:
+    assert(ratio > 1.0)
+
+    downsampled_image_list = list()
+    for image in images:
+        new_spacing = [spacing * ratio for spacing in itk.spacing(image)]
+        new_size = [int(size / ratio) for size in itk.size(image)]
+        downsampled_image = itk.resample_image_filter(image,
+                                                      size=new_size,
+                                                      output_origin=itk.origin(image),
+                                                      output_spacing=new_spacing)
+        downsampled_image_list.append(downsampled_image)
+    return downsampled_image_list
+
+
+# Generate meshes from binary images
+def binary_image_list_to_meshes(images:list, mesh_pixel_type:type=itk.F, object_pixel_value=1) -> list:
+    _, dimension = itk.template(images[0])[1]
+    mesh_type = itk.Mesh[mesh_pixel_type, dimension]
+
+    mesh_list = list()
+    for image in images:
+        mesh = itk.binary_mask3_d_mesh_source(input=image,
+                                              object_value=object_pixel_value,
+                                              ttype=[type(image), mesh_type])
+        mesh_list.append(mesh)
+    return mesh_list
+
+
+# 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):
+    from .pointsetentropyregistrar import PointSetEntropyRegistrar
+
+    registrar = PointSetEntropyRegistrar()
+    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)
+
+    # Run registration and resample from target
+    (transform, deformed_mesh) = registrar.register(template_mesh=template_copy,
+                                                    target_mesh=sample_mesh,
+                                                    metric=metric,
+                                                    transform=transform,
+                                                    learning_rate=learning_rate,
+                                                    max_iterations=max_iterations,
+                                                    verbose=verbose)
+    return deformed_mesh
+
+
+# Align correspondence meshes with Procrustes
+def get_mean_correspondence_mesh(template_meshes:list(),
+                                 convergence_threshold:float=1.0,
+                                 verbose=False) -> itk.Mesh:
+    # Verify ITKShape dependency
+    assert(hasattr(itk, 'MeshProcrustesAlignFilter'))
+
+    # Verify templates have correspondence points
+    assert(all(mesh.GetNumberOfPoints() == template_meshes[0].GetNumberOfPoints()
+               for mesh in template_meshes))
+
+    mesh_type = type(template_meshes[0])
+
+    alignment_filter = itk.MeshProcrustesAlignFilter[mesh_type, mesh_type].New(
+        convergence=convergence_threshold)
+    alignment_filter.SetUseInitialAverageOff()
+    alignment_filter.SetUseNormalizationOff()
+    alignment_filter.SetUseScalingOff()
+    alignment_filter.SetNumberOfInputs(len(template_meshes))
+    for idx in range(len(template_meshes)):
+        alignment_filter.SetInput(idx, template_meshes[idx])
+    alignment_filter.Update()
+
+    if(verbose):
+        print(f'Alignment converged at {alignment_filter.GetMeanPointsDifference()}')
+
+    # Mean output contains only vertices so add cells back before returning
+    mean_mesh_result = alignment_filter.GetMean()
+    mean_mesh_result.SetCells(template_meshes[0].GetCells())
+    return mean_mesh_result
+
+
+# Compute largest correspondence distance between two meshes
+def get_pairwise_hausdorff_distance(first_mesh:itk.Mesh, second_mesh:itk.Mesh) -> float:
+    def euclidean_distance(pt1:itk.Point, pt2:itk.Point) -> float:
+        _, dimension = itk.template(pt1)[1]
+        return sum([(pt1[dim] - pt2[dim]) ** 2 for dim in range(dimension)]) ** 0.5
+
+    max_distance = euclidean_distance(first_mesh.GetPoint(0), second_mesh.GetPoint(0))
+    for idx in range(1, first_mesh.GetNumberOfPoints()):
+        max_distance = max(max_distance,
+                           euclidean_distance(first_mesh.GetPoint(idx), second_mesh.GetPoint(idx)))
+    return max_distance
+
+
+# Refine a template by registering it to a target population,
+# setting template points to nearest neighbors, and then running Procrustes
+# alignment to get the mean of the registered meshes
+def refine_template_from_population(template_mesh:itk.Mesh,
+                                    target_meshes:list,
+                                    registration_iterations=500,
+                                    alignment_threshold=0.1,
+                                    verbose:bool=False) -> itk.Mesh:
+    # Verify ITKShape dependency
+    assert(hasattr(itk, 'MeshProcrustesAlignFilter'))
+
+    # Register template to each target mesh in the population
+    registered_templates = list()
+    for target_mesh in target_meshes:
+        registered_template = register_template_to_sample(template_mesh,
+                                                          target_mesh,
+                                                          max_iterations=registration_iterations,
+                                                          verbose=verbose)
+        registered_templates.append(registered_template)
+
+    # Deform each registered template to correspond with its respective target
+    deformed_templates = list()
+    for idx in range(len(target_meshes)):
+        deformed_templates.append(resample_template_from_target(registered_templates[idx], target_meshes[idx]))
+
+    # Align templates
+    mesh_result = get_mean_correspondence_mesh(deformed_templates,
+                                               convergence_threshold=alignment_threshold,
+                                               verbose=verbose)
+
+    # Quantify difference from original template
+    if(verbose):
+        distance = get_pairwise_hausdorff_distance(mesh_result, template_mesh)
+        print(f'Got Hausdorff distance of {distance} between refined template and original.')
+
+    return mesh_result
diff --git a/src/hasi/hasi/diffeoregistrar.py b/src/hasi/hasi/diffeoregistrar.py
index e3522bd..259c663 100644
--- a/src/hasi/hasi/diffeoregistrar.py
+++ b/src/hasi/hasi/diffeoregistrar.py
@@ -1,7 +1,19 @@
 #!/usr/bin/env python3
 
-# Filename: DiffeoRegistrar.py
-# Contributors: Tom Birdsong
+# Copyright NumFOCUS
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#        http://www.apache.org/licenses/LICENSE-2.0.txt
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
 # Purpose: Python implementation of Slicer SALT mesh-to-mesh registration Cxx
 # CLI.
 #    Ported from https://github.com/slicersalt/RegistrationBasedCorrespondence/
diff --git a/src/hasi/hasi/meansquaresregistrar.py b/src/hasi/hasi/meansquaresregistrar.py
index 1bccc1c..748ddde 100644
--- a/src/hasi/hasi/meansquaresregistrar.py
+++ b/src/hasi/hasi/meansquaresregistrar.py
@@ -1,7 +1,19 @@
 #!/usr/bin/env python3
 
-# Filename: MeanSquaresRegistrar.py
-# Contributors: Tom Birdsong
+# Copyright NumFOCUS
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#        http://www.apache.org/licenses/LICENSE-2.0.txt
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
 # Purpose: Python implementation of Slicer SALT mesh-to-mesh registration Cxx
 # CLI.
 #    Ported from https://github.com/slicersalt/RegistrationBasedCorrespondence/
diff --git a/src/hasi/hasi/meshtomeshregistrar.py b/src/hasi/hasi/meshtomeshregistrar.py
index c56355b..69b044f 100644
--- a/src/hasi/hasi/meshtomeshregistrar.py
+++ b/src/hasi/hasi/meshtomeshregistrar.py
@@ -1,7 +1,19 @@
 #!/usr/bin/env python3
 
-# Filename: MeshToMeshRegistration
-# Contributors: Tom Birdsong
+# Copyright NumFOCUS
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#        http://www.apache.org/licenses/LICENSE-2.0.txt
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
 # Purpose: Base Python class for implementations of mesh-to-mesh registration algorithms.
 #   Initially ported from https://github.com/slicersalt/RegistrationBasedCorrespondence/
 import logging
diff --git a/src/hasi/hasi/pointsetentropyregistrar.py b/src/hasi/hasi/pointsetentropyregistrar.py
index 2d33f3c..f6c3462 100644
--- a/src/hasi/hasi/pointsetentropyregistrar.py
+++ b/src/hasi/hasi/pointsetentropyregistrar.py
@@ -1,11 +1,22 @@
 #!/usr/bin/env python3
 
-# Filename:     pointsetentropyregistrar.py
-# Contributors: Tom Birdsong
+# Copyright NumFOCUS
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#        http://www.apache.org/licenses/LICENSE-2.0.txt
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
 # Purpose:  Python implementation of mesh-to-mesh registration
 #     using point set entropy-based registration
-# Referenced https://github.com/InsightSoftwareConsortium/ITKTrimmedPointSetRegistration/blob/master/examples/TrimmedPointSetRegistrationExample.cxx
-#     for development purposes
+# Related: https://github.com/InsightSoftwareConsortium/ITKTrimmedPointSetRegistration/blob/master/examples/TrimmedPointSetRegistrationExample.cxx
 
 import logging
 import os


Test patch:
diff --git a/.github/workflows/test-python-hasi.yml b/.github/workflows/test-python-hasi.yml
index fac5da5..4e7ee79 100644
--- a/.github/workflows/test-python-hasi.yml
+++ b/.github/workflows/test-python-hasi.yml
@@ -28,8 +28,8 @@ jobs:
     - name: Install build dependencies
       run: |
         python -m pip install --upgrade pip
-        python -m pip install ninja
         python -m pip install itk==5.2.0
+        python -m pip install itk-shape==0.2.0
         python -m pip install pytest
 
     - name: Test with pytest
diff --git a/src/hasi/test_align.py b/src/hasi/test_align.py
new file mode 100644
index 0000000..f2f789a
--- /dev/null
+++ b/src/hasi/test_align.py
@@ -0,0 +1,185 @@
+#!/usr/bin/env python3
+
+# Copyright NumFOCUS
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#        http://www.apache.org/licenses/LICENSE-2.0.txt
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+# Purpose:  Pytest definitions for atlas alignment classes
+#           defined in the src/hasi/hasi module
+
+# Directions:
+# - run `pytest .`
+
+import os
+import sys
+from urllib.request import urlretrieve
+module_path = os.path.abspath(os.path.join('.'))
+
+if module_path not in sys.path:
+  sys.path.append(module_path)
+
+import itk
+
+MEANSQUARES_METRIC_MAXIMUM_THRESHOLD = 0.015
+DIFFEO_METRIC_MAXIMUM_THRESHOLD = 0.25
+POINT_SET_METRIC_MAXIMUM_THRESHOLD = 0.0001
+MAX_ITERATIONS = 200
+
+TEMPLATE_MESH_FILE = 'test/Input/906-R-atlas.obj'
+TARGET_MESH_FILE = 'test/Input/901-R-mesh.vtk'
+TARGET_HEALTHY_IMAGE_FILE = 'test/Input/901-R-femur-label.nrrd'
+TARGET_UNHEALTHY_IMAGE_FILE = 'test/Input/901-L-femur-label.nrrd'
+
+os.makedirs('test', exist_ok=True)
+os.makedirs('test/Input', exist_ok=True)
+os.makedirs('test/Output', exist_ok=True)
+
+if not os.path.exists(TEMPLATE_MESH_FILE):
+    url = 'https://data.kitware.com/api/v1/file/608b006d2fa25629b970f139/download'
+    urlretrieve(url, TEMPLATE_MESH_FILE)
+    
+if not os.path.exists(TARGET_MESH_FILE):
+    url = 'https://data.kitware.com/api/v1/file/5f9daaba50a41e3d1924dae9/download'
+    urlretrieve(url, TARGET_MESH_FILE)
+    
+if not os.path.exists(TARGET_HEALTHY_IMAGE_FILE):
+    url = 'https://data.kitware.com/api/v1/file/5eea20bd9014a6d84ec75df1/download'
+    urlretrieve(url, TARGET_HEALTHY_IMAGE_FILE)
+    
+if not os.path.exists(TARGET_UNHEALTHY_IMAGE_FILE):
+    url = 'https://data.kitware.com/api/v1/file/5eea20bb9014a6d84ec75de5/download'
+    urlretrieve(url, TARGET_UNHEALTHY_IMAGE_FILE)
+    
+template_mesh = itk.meshread(TEMPLATE_MESH_FILE, itk.F)
+target_mesh = itk.meshread(TARGET_MESH_FILE, itk.F)
+target_healthy_image = itk.imread(TARGET_HEALTHY_IMAGE_FILE)
+target_unhealthy_image = itk.imread(TARGET_UNHEALTHY_IMAGE_FILE)
+
+mesh_pixel_type, dimension = itk.template(template_mesh)[1]
+
+def test_hasi_import():
+    # Import succeeds
+    from hasi import align
+
+
+def test_mesh_to_image():
+    # Import succeeds
+    from hasi.align import mesh_to_image
+
+    # Mesh-to-image executes without error
+    images = mesh_to_image([template_mesh, target_mesh])
+
+    # Images were generated
+    assert(all([type(image) == itk.Image[mesh_pixel_type, dimension] for image in images]))
+
+    # Images occupy the same physical space
+    assert(all([itk.spacing(image) == itk.spacing(images[0]) for image in images]))
+    assert(all([itk.size(image) == itk.size(images[0]) for image in images]))
+    assert(all([itk.origin(image) == itk.origin(images[0]) for image in images]))
+
+    # Mesh-to-image can run with reference image
+    image_from_reference = mesh_to_image(target_mesh, reference_image=images[0])
+
+    # Type and image attributes match previous output
+    assert(type(image_from_reference) == type(images[0]))
+    assert(itk.spacing(image_from_reference) == itk.spacing(images[0]))
+    assert(itk.size(image_from_reference) == itk.size(images[0]))
+    assert(itk.origin(image_from_reference) == itk.origin(images[0]))
+
+
+def test_resample_from_target():
+    DISTANCE_THRESHOLD = 0.7
+
+    # Import succeeds
+    from hasi.align import resample_template_from_target, get_pairwise_hausdorff_distance
+
+    # Resample succeeds
+    # Note for actual use we would want to register first
+    resampled_template = resample_template_from_target(template_mesh, target_mesh)
+
+    # Template has same number of points before and after
+    assert(resampled_template.GetNumberOfPoints() == template_mesh.GetNumberOfPoints())
+
+    # Template points moved
+    assert(resampled_template.GetPoint(0) != template_mesh.GetPoint(0))
+
+    # All template points lie on template mesh
+    assert(any([resampled_template.GetPoint(1) == target_mesh.GetPoint(j)
+                    for j in range(target_mesh.GetNumberOfPoints())]))
+
+    # Movement distance was within threshold
+    distance = get_pairwise_hausdorff_distance(resampled_template, template_mesh)
+    assert(distance > 0.0 and distance < DISTANCE_THRESHOLD)
+
+
+def test_paste_to_common_space():
+    # Import succeeds
+    from hasi.align import mesh_to_image, paste_to_common_space
+    
+    # Paste succeeds
+    resized_images = paste_to_common_space([target_healthy_image, target_unhealthy_image])
+    assert(len(resized_images) == 2)
+
+    # Verify common space
+    TOLERANCE = 1e-6
+    assert(itk.size(resized_images[0]) == itk.size(resized_images[1]))
+    assert(all([itk.spacing(resized_images[0])[dim] - itk.spacing(resized_images[1])[dim] < TOLERANCE
+                for dim in range(dimension)]))
+
+
+def test_downsample_images():
+    DOWNSAMPLE_RATIO = 2.0
+
+    # Import succeeds
+    from hasi.align import mesh_to_image, downsample_images
+
+    # Downsample succeeds
+    images = downsample_images([target_healthy_image, target_unhealthy_image], DOWNSAMPLE_RATIO)
+    assert(len(images) == 2)
+
+    # Downsampled image dimensions adjusted by expected ratio
+    assert(itk.size(images[0]) ==
+           [int(itk.size(target_healthy_image)[dim] / DOWNSAMPLE_RATIO) for dim in range(dimension)])
+    assert(itk.spacing(images[0]) ==
+           [itk.spacing(target_healthy_image)[dim] * DOWNSAMPLE_RATIO for dim in range(dimension)])
+    assert(itk.size(images[1]) ==
+           [int(itk.size(target_unhealthy_image)[dim] / DOWNSAMPLE_RATIO) for dim in range(dimension)])
+    assert(itk.spacing(images[1]) ==
+           [itk.spacing(target_unhealthy_image)[dim] * DOWNSAMPLE_RATIO for dim in range(dimension)])
+
+
+def test_generate_meshes():
+    # Import succeeds
+    from hasi.align import binary_image_list_to_meshes
+
+    # Mesh generation succeeds
+    meshes = binary_image_list_to_meshes([target_healthy_image, target_unhealthy_image], itk.F, 1)
+    assert(len(meshes) == 2)
+
+    # Mesh vertices match expectation
+    assert(meshes[0].GetNumberOfPoints() == 900671)
+    assert(meshes[1].GetNumberOfPoints() == 855276)
+
+
+def test_refine_template_from_population():
+    # Import succeeds
+    from hasi.align import refine_template_from_population, get_pairwise_hausdorff_distance
+
+    # Refinement succeeds
+    mesh_result = refine_template_from_population(template_mesh,
+                                                  [target_mesh],
+                                                  registration_iterations=100)
+
+    # Distance matches expectation
+    distance = get_pairwise_hausdorff_distance(mesh_result, template_mesh)
+    assert(distance > 0.7 and distance < 1.0)
diff --git a/src/hasi/test_meshtomeshregistration.py b/src/hasi/test_registrars.py
similarity index 87%
rename from src/hasi/test_meshtomeshregistration.py
rename to src/hasi/test_registrars.py
index 151880c..9ad476f 100644
--- a/src/hasi/test_meshtomeshregistration.py
+++ b/src/hasi/test_registrars.py
@@ -1,12 +1,23 @@
 #!/usr/bin/env python3
 
-# Filename:     test_meshtomeshregistration.py
-# Contributors: Tom Birdsong
+# Copyright NumFOCUS
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#        http://www.apache.org/licenses/LICENSE-2.0.txt
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
 # Purpose:  Pytest definitions for mesh-to-mesh registration classes
 #           defined in the src/hasi/hasi module
 
 # Directions:
-# - cd into the src/hasi directory
 # - run `pytest .`
 
 import os
@@ -42,31 +53,6 @@
 template_mesh = itk.meshread(TEMPLATE_MESH_FILE, itk.F)
 target_mesh = itk.meshread(TARGET_MESH_FILE, itk.F)
 
-# Test base class import and functions
-def test_mesh_to_image():
-    # Class is imported
-    from hasi.meshtomeshregistrar import MeshToMeshRegistrar
-
-    # Class is instantiable
-    registrar = MeshToMeshRegistrar()
-
-    # Initialize runs without error
-    registrar.initialize()
-
-    # Mesh-to-image executes without error
-    img1 = registrar.mesh_to_image(template_mesh)
-
-    # Image was generated
-    assert(type(img1) == itk.Image[itk.F, 3])
-
-    # Mesh-to-image can run with reference image
-    img2 = registrar.mesh_to_image(target_mesh, img1)
-    assert(type(img2) == itk.Image[itk.F, 3])
-
-    # Images occupy the same physical space
-    assert(img2.GetOrigin() == img1.GetOrigin())
-    assert(img2.GetSpacing() == img1.GetSpacing())
-
 # Test registration with mean squares image metric
 def test_meansquares_registration():
     MESH_OUTPUT_FILE = 'test/Output/testMeanSquaresRegisterOutput.vtk'