> ## Documentation Index
> Fetch the complete documentation index at: https://docs.ntop.com/llms.txt
> Use this file to discover all available pages before exploring further.

# Mesh (Python)

This example loads a `.implicit` file (a heat sink), samples the body's field over a regular 3D grid using the batched `ntop_core_query_field_array` call, and reconstructs a surface mesh from the sampled values using scikit-image's `measure.marching_cubes`. The resulting mesh is written to an STL file and also rendered to a PNG with Matplotlib.

```python theme={null}
"""A simple implementation of the "Marching Cubes" using scikit-image."""
# region Imports
from ctypes import (
    CDLL,
    CFUNCTYPE,
    pointer,
    POINTER,
    c_char_p,
    c_double,
    c_uint32,
    c_void_p,
    Structure
)
import os
from pathlib import Path
from typing import Any

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import numpy as np
import numpy.typing as npt
from skimage import measure
# endregion Imports

#########
#########
# Setup #
#########
#########
# region Preamble Setup

########################
# nTop Core Primitives #
########################
# region nTop Core Primitives

class Vec3(Structure):
    """X, Y, Z components representing a point or vector."""
    _fields_ = [("x", c_double), ("y", c_double), ("z", c_double)]

class BoundingBox(Structure):
    """
    Minimal, axis-aligned bounding box of an Implicit.

    The representation is comprised of two points (min_x, min_y, min_z) and (max_x, max_y, max_z).
    """
    _fields_ = [("min", Vec3), ("max", Vec3)]

# endregion

###########################################################
# Setup library functions and data marshaling definitions #
###########################################################
# region Library wrapping setup 

# Load the nTop Core library
NTOP_CORE_LIB_ENVVAR = "NTOP_CORE_LIB"
if NTOP_CORE_LIB_ENVVAR in os.environ:
    env_lib_path = Path(os.environ[NTOP_CORE_LIB_ENVVAR])
    if env_lib_path.exists() and env_lib_path.is_file():
        path_to_lib = env_lib_path
    else:
        raise FileNotFoundError(env_lib_path)

ntop_core = CDLL(str(path_to_lib))
# In order to successfully load the nTop Core Library the dependent libraries of
# the `ntop_core` library must be on the PATH or the current working directory must contain
# those libraries not found on the PATH.
# For simplicity here we are choosing, for the moment, to append to the PATH environment
if str(path_to_lib.parent) not in os.environ["PATH"]:
    os.environ["PATH"] += os.pathsep + str(path_to_lib.parent)

# Define data marshalling specification for loading implicit
ntop_core.ntop_core_import_from_file.restype = c_uint32
ntop_core.ntop_core_import_from_file.argtypes = [
    c_char_p,
    POINTER(c_void_p)
]

def callback_import_from_file_error_check(result: c_uint32, func, arguments: Any) -> bool:
    """
    Callback function which happens after the function is called but before 
    control is returned to the caller enabling error handling.

    :param result: Return value of `ntop_core_import_from_file`
    :param arguments: Tuple of parameters passed to the invocation
        of `ntop_core_import_from_file`
    :returns: True if loading the implicit was successful.
    :raises: RuntimeError if any error state is reported by nTop Core.
    """
    # Check to see if the implicit loaded correctly
    if result != 0:
        raise RuntimeError(f"Cannot load implicit: {arguments[0]}.")

    return True

# Setup error check when loading the file.
ntop_core.ntop_core_import_from_file.errcheck = callback_import_from_file_error_check

# Define data marshalling specification for releasing the memory used for the `.implicit`
ntop_core.ntop_core_release.restype = None
ntop_core.ntop_core_release.argtypes = [c_void_p]

# Define data marshalling specification for getting the bounding box of an implicit
ntop_core.ntop_core_query_bounding_box.restype = None
ntop_core.ntop_core_query_bounding_box.argtypes = [
    c_void_p,
    POINTER(BoundingBox)
]

# Define the signature for the callback in the `ntop_core_query_field_array``
CFUNCTYPE_QUERY_FIELD_ARRAY_CALLBACK = CFUNCTYPE(
    None,
    c_void_p,
    POINTER(c_double),
    c_uint32
)

# An example callback that can be used with `ntop_core_query_field_array` invocations.
def query_field_array_callback(scope_pass) -> CFUNCTYPE_QUERY_FIELD_ARRAY_CALLBACK:
    """Callback for the ntop_core_query_field_array to smuggle out the query data."""

    def _smuggle(context: c_void_p, query_vals: POINTER(c_double), count: int) -> None:
        temp = np.ctypeslib.as_array(query_vals, shape = (count, 1)).flatten()
        np.copyto(scope_pass, temp)

    return  CFUNCTYPE_QUERY_FIELD_ARRAY_CALLBACK(_smuggle)

ntop_core.ntop_core_query_field_array.restype = None
ntop_core.ntop_core_query_field_array.argtypes = [
    c_void_p,
    np.ctypeslib.ndpointer(dtype=np.double, ndim=2, flags='CONTIGUOUS'),
    c_uint32,
    c_void_p,
    CFUNCTYPE_QUERY_FIELD_ARRAY_CALLBACK
]
# endregion Library wrapping setup 

# endregion Preamble Setup

######################
######################
# Utility Functions #
######################
######################
# region Utilities

###################
# STL File Writer #
###################
# region STL writer
def write_stl(triangles: npt.ArrayLike, filename: Path) -> None:
    """
    Write a mesh stored as triangles in STL format.
    
    Note: degenerate faces where vertices are colinear are filtered out.

    :param triangles: numpy Array of dimensions Nx3x3 (N is the number of triangles in
                    the mesh) with each vertex in 3D space. 
    :param filename: File name of the STL to write.
    """
    with open(filename, "w") as f:
        f.write("solid object\n")
        for triangle in triangles:
            normal = np.cross(triangle[1] - triangle[0], triangle[2] - triangle[0])
            normal_length = np.sqrt(np.sum(normal ** 2))
            if (normal_length == 0):
                # Skip degenerate faces
                continue

            unit_normal = normal / normal_length

            f.write(f"  facet normal {unit_normal[0]:.6f}  {unit_normal[1]:.6f}  {unit_normal[2]:.6f}\n")
            f.write(f"    outer loop\n")
            for point in triangle:
                f.write(f"      vertex {point[0]:.6f}  {point[1]:.6f}  {point[2]:.6f}\n")
            f.write(f"    endloop\n")
            f.write(f"  endfacet\n")
        f.write("endsolid object\n")

# endregion STL writer

# endregion Utilities

################################################
################################################
# Functional example of meshing with nTop Core #
################################################
#################################################
# region Example meshing code with nTop Core

# Load the implicit
gyroid_implicit_file_path = Path(__file__).parent.parent / "assets" / "heat-sink.implicit"
gyroid_handle = c_void_p()
ntop_core.ntop_core_import_from_file(str(gyroid_implicit_file_path).encode("utf8"), gyroid_handle)

# Get the bounding box to know the extent of the points to mesh
bounding_box = BoundingBox()
ntop_core.ntop_core_query_bounding_box(gyroid_handle, pointer(bounding_box))

# Offset the bounding box slightly to ensure that we capture points beyond the isosurface
# so the meshing algorithm can close the part.
bounding_box_offset = 0.001

# Set the minimum feature size you want to capture in the part. The values
# selected here will have an O(n^3) impact on mesh file size amd generation
# time as it directly effects the number of sampled points.
# In this example it is isotropic but if dimensional variability is important
# that can be modified here.
min_feature_size = (0.00025, 0.00025, 0.00025)

# Setup the sampling grid
sample_points_numpy_form = np.mgrid[
    bounding_box.min.x - bounding_box_offset:bounding_box.max.x + bounding_box_offset:min_feature_size[0],
    bounding_box.min.y - bounding_box_offset:bounding_box.max.y + bounding_box_offset:min_feature_size[1],
    bounding_box.min.z - bounding_box_offset:bounding_box.max.z + bounding_box_offset:min_feature_size[2],
]

# Turn this into an array of x,y,z points as expected by nTop Core
sample_points_ntop_core_form = np.ascontiguousarray(sample_points_numpy_form.reshape(3, -1, order='C').transpose(), dtype=np.double)

# Setup the variable to capture the output from nTop Core
field_values = np.zeros(sample_points_ntop_core_form.shape[0])

# Use nTop Core's batch query functionality to query the field at all the
# points simultaneously.
# Note: for sufficiently large parts or very fine features you may need to
# chunk these queries depending on your system configuration.
ntop_core.ntop_core_query_field_array(
    gyroid_handle,
    sample_points_ntop_core_form,
    sample_points_ntop_core_form.shape[0],
    None,
    query_field_array_callback(field_values)
)

# Free the memory associated with the implicit as no more calls are needed
ntop_core.ntop_core_release(gyroid_handle)

# Arrange the data back to the original shape as that is the form the meshing
# algorithm expects: a top-down Z-axis slice view of field values.
volume_data = field_values.reshape(sample_points_numpy_form.shape[1:])

# Use skimage to mesh the data
vertices, faces, normals, values = measure.marching_cubes(
    volume_data,
    0,
    allow_degenerate=False,
    spacing=min_feature_size
)

# The implementation of `measure.marching_cubes()` translates samples to its
# own frame in octant 1 so we translate to the original position of the implicit
# for consistency.
vertices += np.array([bounding_box.min.x, bounding_box.min.y, bounding_box.min.z])
vertices -= bounding_box_offset

# Fancy indexing: `vertices[faces]` to generate a collection of triangles.
# This takes the form of an array containing triplets of 3D points.
triangles = vertices[faces]

# Make a folder to capture the output
output_folder = Path(__file__).parent / "output"
output_folder.mkdir(exist_ok=True)

# Write the mesh as an STL
generated_mesh_filepath = output_folder / "mesh.stl"
write_stl(triangles, generated_mesh_filepath)

# Display resulting mesh using Matplotlib.
fig = plt.figure(figsize=(25, 25))
ax = fig.add_subplot(111, projection='3d')

# Add the mesh visualization
# Note: Matplotlib is not very performant for this rendering task and should
# primarily be used for proof-of-concept using other mesh visualization tools
# for deeper analysis.
mesh = Poly3DCollection(triangles)
mesh.set_edgecolor('k')
ax.add_collection3d(mesh)

# Label the plot
ax.set_title("Mesh visualization of implicit")
ax.set_xlabel("x-axis (meters)")
ax.set_ylabel("y-axis (meters)")
ax.set_zlabel("z-axis (meters)")

# Set reasonable limits to scale the part appropriately for viewing
ax.set_xlim(bounding_box.min.x, bounding_box.max.x)
ax.set_ylim(bounding_box.min.y, bounding_box.max.y)
ax.set_zlim(bounding_box.min.z, bounding_box.max.z)

# Render the mesh and save to a PNG
mesh_render_filepath = output_folder / "mesh.png"
plt.savefig(mesh_render_filepath)
print(f"Render of mesh saved to {mesh_render_filepath}")

# Uncomment the line below to render the mesh in an interactive window
# plt.show()

# endregion Example meshing code with nTop Core
```