diff --git a/itk-module.cmake b/itk-module.cmake
index 1c1b3c1..d1dd663 100644
--- a/itk-module.cmake
+++ b/itk-module.cmake
@@ -27,6 +27,10 @@ itk_module(Cuberille
     ITKQuadEdgeMeshFiltering
     ITKThresholding
     ITKIOImageBase
+    ITKIOMeta
+    ITKIONRRD
+    ITKIOMeshBase
+    ITKIOVTK
   DESCRIPTION
     "${DOCUMENTATION}"
   EXCLUDE_FROM_DEFAULT
diff --git a/wrapping/test/CMakeLists.txt b/wrapping/test/CMakeLists.txt
new file mode 100644
index 0000000..026381e
--- /dev/null
+++ b/wrapping/test/CMakeLists.txt
@@ -0,0 +1,41 @@
+itk_python_expression_add_test(
+  NAME itkCuberilleImageToMeshFilterPythonTest
+  EXPRESSION "filter = itk.CuberilleImageToMeshFilter.New()"
+)
+
+set(test_input_dir ${CMAKE_CURRENT_SOURCE_DIR}/../../test/Input)
+# cmake_path(CONVERT ${test_input_dir} TO_CMAKE_PATH_LIST test_input_dir NORMALIZE)  # CMake 3.20+
+
+list(FIND ITK_WRAP_IMAGE_DIMS 3 wrap_3_index)
+if(wrap_3_index GREATER -1)
+  itk_python_add_test(NAME Cuberille_blob0_python
+    TEST_DRIVER_ARGS
+    COMMAND cuberille_test.py
+      DATA{${test_input_dir}/blob0.mha}
+      ${ITK_TEST_OUTPUT_DIR}/blob0.vtk
+    )
+
+  itk_python_add_test(NAME Cuberille_blob0_100_python
+    TEST_DRIVER_ARGS
+    COMMAND cuberille_test.py
+      DATA{${test_input_dir}/blob0.mha}
+      ${ITK_TEST_OUTPUT_DIR}/blob0_100.vtk
+        100  # iso-surface value
+    )
+
+  itk_python_add_test(NAME Cuberille_hydrogenAtom_35_python
+    TEST_DRIVER_ARGS
+    COMMAND cuberille_test.py
+      DATA{${test_input_dir}/hydrogenAtom.mha}
+      ${ITK_TEST_OUTPUT_DIR}/hydrogenAtom35.vtk
+        35  # iso-surface value
+    )
+
+  itk_python_add_test(NAME Cuberille_hydrogenAtom_50_python
+    TEST_DRIVER_ARGS
+    COMMAND cuberille_test.py
+      DATA{${test_input_dir}/hydrogenAtom.mha}
+      ${ITK_TEST_OUTPUT_DIR}/hydrogenAtom50.vtk
+        50  # iso-surface value
+    )
+endif()
diff --git a/wrapping/test/cuberille_test.py b/wrapping/test/cuberille_test.py
new file mode 100644
index 0000000..ee90a25
--- /dev/null
+++ b/wrapping/test/cuberille_test.py
@@ -0,0 +1,30 @@
+import sys
+
+import itk
+
+image = itk.imread(sys.argv[1])
+
+if len(sys.argv) > 3:
+    isosurface = sys.argv[3]
+else:
+    calculator = itk.MinimumMaximumImageCalculator[type(image)].New()
+    calculator.SetImage(image)
+    calculator.Compute()
+    isosurface = (calculator.GetMaximum() + calculator.GetMinimum()) / 2
+
+# convert number type
+(input_image_template, (input_pixel_type, input_image_dimension)) = itk.template(image)
+if input_pixel_type == itk.F or input_pixel_type == itk.D:
+    isosurface = float(isosurface)
+else:
+    isosurface = int(isosurface)
+
+mesh = itk.cuberille_image_to_mesh_filter(
+    image,
+    generate_triangle_faces=True,
+    iso_surface_value=isosurface,
+    project_vertices_to_iso_surface=True,
+    project_vertex_surface_distance_threshold=0.05,
+)
+
+itk.meshwrite(mesh, sys.argv[2])