Mesh Creation: Building, Smoothing, and Refining the Brain Mesh¶

Overview¶

This notebook takes the labelled voxel-based geometry produced by 02_ImageProcess.ipynb and converts it into a finite element mesh.

The main steps are:

0. set up the subject paths, mesh parameters, and mesh smoother executable;
1. convert pre_model.nii.gz into a hexahedral mesh and add model structures such as the falx, tentorium, pia, and dura;
2. smooth the mesh surfaces to reduce voxel stair-step artefacts;
3. repair brain-skull contact so that the brain is separated from the skull by CSF.

The expected input is:

data/subjects/<subject_name>/pre_model.nii.gz

The final outputs are saved in:

data/subjects/<subject_name>/output/

0. Setup paths and parameters¶

This section defines the subject, project folders, dependency paths, and mesh settings used throughout the notebook.

The most important input is pre_model.nii.gz, which is the labelled geometry created in the previous image-processing notebook. The main output folder is output/ inside the subject folder. This is where the raw mesh, smoothed mesh, and revised mesh will be written.

In this step, users should check the following settings before running the notebook:

• subject_name: the subject folder to process.
• mesh_size: the target voxel/element size used for mesh generation.
• iterations: the maximum number of local refinement iterations during mesh creation.
• blocks: the number of processing blocks used when generating the mesh.
• threshold: the threshold used when converting the labelled image into the mesh.
• smoother_name: the mesh smoother executable used for smoothing the generated mesh.
• smoothing parameters such as smoothing_passes, relaxation_factor, surface_iteration_limit, and internal_search_limit.

The mesh smoother executable must match the user's operating system and processor type. For example, a macOS Apple Silicon executable will not run on Windows or Linux, and a Linux executable will not run directly on macOS. Users should select the correct executable name in Step 0 before running the smoothing step.

For most users, the most important values to check are subject_name, mesh_size, and smoother_name. A smaller mesh_size produces a finer mesh but increases the number of elements and computational cost.

In [2]:

from pathlib import Path
import os
import sys

# --------------------------------------------------
# User settings
# --------------------------------------------------
# Change this value when running the notebook for a different subject.
subject_name = "sub0045"

# Name of the folder created inside the subject directory for mesh outputs.
# The mesh creation code writes files such as mesh.k and mesh_smoothed.k here.
output_name = "output"

# Target mesh size in mm. Smaller values preserve more detail but generate
# larger meshes; larger values are faster but may reduce geometric detail.
mesh_size = 2

# --------------------------------------------------
# Voxel correction settings
# --------------------------------------------------
# These parameters control the custom voxel-cleaning step before mesh creation.
# In most cases, they should be left unchanged.
iterations = 3      # Number of correction passes through the labelled geometry.
blocks = 1          # Neighbourhood size used when checking surrounding voxels.
threshold = 0.9     # Neighbour agreement threshold for relabelling uncertain voxels.

# --------------------------------------------------
# Optional centre-of-gravity settings
# --------------------------------------------------
# These values are kept for downstream simulation setup. They are not directly
# used by the standard mesh creation commands below unless explicitly added to
# later processing steps.
cog_id = 80000000
change_ap = 0       # Anterior-posterior shift.
change_lr = 0       # Left-right shift.
change_is = 0       # Inferior-superior shift.
change_cog = [change_ap, change_lr, change_is]

# --------------------------------------------------
# Project paths
# --------------------------------------------------
# This assumes the notebook is run from the repository's notebooks/ folder.
PROJECT_ROOT = Path.cwd().resolve().parent

# Source code folder containing the brain_mesh_creation Python package.
SRC_DIR = PROJECT_ROOT / "src"
PACKAGE_DIR = SRC_DIR / "brain_mesh_creation"

# Dependency folder containing reference files, part lists, and executables.
DEPENDENCIES_DIR = SRC_DIR / "dependencies"
RS_DIR = DEPENDENCIES_DIR / "rs"

# Subject data folder.
SUBJECTS_DIR = PROJECT_ROOT / "data" / "subjects"
SUBJECT_DIR = SUBJECTS_DIR / subject_name

# Input geometry from 02_ImageProcess.ipynb.
PRE_MODEL_FILE = SUBJECT_DIR / "pre_model.nii.gz"

# Mesh output folder created by the mesh generation step.
OUTPUT_DIR = SUBJECT_DIR / output_name

# Mesh smoother executable. Change this filename if running on another platform.
# Common options in src/dependencies/rs/ are expected to include:
#   mesh_smoother_linux_x86_64
#   mesh_smoother_macos_x86_64
#   mesh_smoother_macos_arm64
#   mesh_smoother_windows_x86_64.exe
smoother_executable = "mesh_smoother_macos_arm64"
SMOOTHER_PATH = RS_DIR / smoother_executable

# --------------------------------------------------
# Make the local package importable
# --------------------------------------------------
# This allows the notebook to import brain_mesh_creation without installing it
# globally into the Python environment.
if str(SRC_DIR) not in sys.path:
    sys.path.insert(0, str(SRC_DIR))

from brain_mesh_creation import bmctk
from brain_mesh_creation import mesh_utils

# --------------------------------------------------
# Export paths for bash cells
# --------------------------------------------------
# Bash cells cannot directly see Python pathlib objects, so the key paths and
# settings are exported as environment variables.
env_vars = {
    "PROJECT_ROOT": PROJECT_ROOT,
    "SRC_DIR": SRC_DIR,
    "PACKAGE_DIR": PACKAGE_DIR,
    "DEPENDENCIES_DIR": DEPENDENCIES_DIR,
    "RS_DIR": RS_DIR,
    "SUBJECTS_DIR": SUBJECTS_DIR,
    "SUBJECT_DIR": SUBJECT_DIR,
    "PRE_MODEL_FILE": PRE_MODEL_FILE,
    "OUTPUT_DIR": OUTPUT_DIR,
    "SMOOTHER_PATH": SMOOTHER_PATH,
}

for key, value in env_vars.items():
    os.environ[key] = str(value)

os.environ["SUBJECT_NAME"] = subject_name
os.environ["OUTPUT_NAME"] = output_name
os.environ["MESH_SIZE"] = str(mesh_size)

# --------------------------------------------------
# Setup summary
# --------------------------------------------------
print("Setup complete.")
print(f"Subject: {subject_name}")
print(f"Mesh size set to: {mesh_size} mm")

from pathlib import Path import os import sys # -------------------------------------------------- # User settings # -------------------------------------------------- # Change this value when running the notebook for a different subject. subject_name = "sub0045" # Name of the folder created inside the subject directory for mesh outputs. # The mesh creation code writes files such as mesh.k and mesh_smoothed.k here. output_name = "output" # Target mesh size in mm. Smaller values preserve more detail but generate # larger meshes; larger values are faster but may reduce geometric detail. mesh_size = 2 # -------------------------------------------------- # Voxel correction settings # -------------------------------------------------- # These parameters control the custom voxel-cleaning step before mesh creation. # In most cases, they should be left unchanged. iterations = 3 # Number of correction passes through the labelled geometry. blocks = 1 # Neighbourhood size used when checking surrounding voxels. threshold = 0.9 # Neighbour agreement threshold for relabelling uncertain voxels. # -------------------------------------------------- # Optional centre-of-gravity settings # -------------------------------------------------- # These values are kept for downstream simulation setup. They are not directly # used by the standard mesh creation commands below unless explicitly added to # later processing steps. cog_id = 80000000 change_ap = 0 # Anterior-posterior shift. change_lr = 0 # Left-right shift. change_is = 0 # Inferior-superior shift. change_cog = [change_ap, change_lr, change_is] # -------------------------------------------------- # Project paths # -------------------------------------------------- # This assumes the notebook is run from the repository's notebooks/ folder. PROJECT_ROOT = Path.cwd().resolve().parent # Source code folder containing the brain_mesh_creation Python package. SRC_DIR = PROJECT_ROOT / "src" PACKAGE_DIR = SRC_DIR / "brain_mesh_creation" # Dependency folder containing reference files, part lists, and executables. DEPENDENCIES_DIR = SRC_DIR / "dependencies" RS_DIR = DEPENDENCIES_DIR / "rs" # Subject data folder. SUBJECTS_DIR = PROJECT_ROOT / "data" / "subjects" SUBJECT_DIR = SUBJECTS_DIR / subject_name # Input geometry from 02_ImageProcess.ipynb. PRE_MODEL_FILE = SUBJECT_DIR / "pre_model.nii.gz" # Mesh output folder created by the mesh generation step. OUTPUT_DIR = SUBJECT_DIR / output_name # Mesh smoother executable. Change this filename if running on another platform. # Common options in src/dependencies/rs/ are expected to include: # mesh_smoother_linux_x86_64 # mesh_smoother_macos_x86_64 # mesh_smoother_macos_arm64 # mesh_smoother_windows_x86_64.exe smoother_executable = "mesh_smoother_macos_arm64" SMOOTHER_PATH = RS_DIR / smoother_executable # -------------------------------------------------- # Make the local package importable # -------------------------------------------------- # This allows the notebook to import brain_mesh_creation without installing it # globally into the Python environment. if str(SRC_DIR) not in sys.path: sys.path.insert(0, str(SRC_DIR)) from brain_mesh_creation import bmctk from brain_mesh_creation import mesh_utils # -------------------------------------------------- # Export paths for bash cells # -------------------------------------------------- # Bash cells cannot directly see Python pathlib objects, so the key paths and # settings are exported as environment variables. env_vars = { "PROJECT_ROOT": PROJECT_ROOT, "SRC_DIR": SRC_DIR, "PACKAGE_DIR": PACKAGE_DIR, "DEPENDENCIES_DIR": DEPENDENCIES_DIR, "RS_DIR": RS_DIR, "SUBJECTS_DIR": SUBJECTS_DIR, "SUBJECT_DIR": SUBJECT_DIR, "PRE_MODEL_FILE": PRE_MODEL_FILE, "OUTPUT_DIR": OUTPUT_DIR, "SMOOTHER_PATH": SMOOTHER_PATH, } for key, value in env_vars.items(): os.environ[key] = str(value) os.environ["SUBJECT_NAME"] = subject_name os.environ["OUTPUT_NAME"] = output_name os.environ["MESH_SIZE"] = str(mesh_size) # -------------------------------------------------- # Setup summary # -------------------------------------------------- print("Setup complete.") print(f"Subject: {subject_name}") print(f"Mesh size set to: {mesh_size} mm")

Setup complete.
Subject: sub0045
Mesh size set to: 2 mm

In [3]:

# --------------------------------------------------
# Check required files and folders
# --------------------------------------------------
required_paths = {
    "Subject folder": SUBJECT_DIR,
    "Input pre-model geometry": PRE_MODEL_FILE,
    "Dependency folder": DEPENDENCIES_DIR,
    "Resource/executable folder": RS_DIR,
    "Mesh smoother executable": SMOOTHER_PATH,
}

missing_paths = [name for name, path in required_paths.items() if not path.exists()]

if missing_paths:
    print("Missing required paths:")
    for name in missing_paths:
        print(f"  - {name}: {required_paths[name]}")
    raise FileNotFoundError("One or more required files or folders are missing.")

print("All required files and folders found.")

# -------------------------------------------------- # Check required files and folders # -------------------------------------------------- required_paths = { "Subject folder": SUBJECT_DIR, "Input pre-model geometry": PRE_MODEL_FILE, "Dependency folder": DEPENDENCIES_DIR, "Resource/executable folder": RS_DIR, "Mesh smoother executable": SMOOTHER_PATH, } missing_paths = [name for name, path in required_paths.items() if not path.exists()] if missing_paths: print("Missing required paths:") for name in missing_paths: print(f" - {name}: {required_paths[name]}") raise FileNotFoundError("One or more required files or folders are missing.") print("All required files and folders found.")

All required files and folders found.

Optional: remove an existing mesh output folder¶

Run the next cell only if you have already generated a mesh for this subject and want to start again from a clean output folder.

By default, the cleanup flag is set to False so the cell will not delete anything accidentally.

In [5]:

import shutil

# Set this to True only when you intentionally want to delete the existing
# output folder and regenerate the mesh from scratch.
RUN_CLEANUP = False

print(f"Subject: {subject_name}")
print(f"Output folder: {OUTPUT_DIR}")

if RUN_CLEANUP:
    if OUTPUT_DIR.exists():
        shutil.rmtree(OUTPUT_DIR)
        print(f"Removed existing output folder. ")
    else:
        print("No existing output folder found. Nothing to remove.")
else:
    print("Cleanup skipped. Set RUN_CLEANUP = True to remove the existing output folder.")

import shutil # Set this to True only when you intentionally want to delete the existing # output folder and regenerate the mesh from scratch. RUN_CLEANUP = False print(f"Subject: {subject_name}") print(f"Output folder: {OUTPUT_DIR}") if RUN_CLEANUP: if OUTPUT_DIR.exists(): shutil.rmtree(OUTPUT_DIR) print(f"Removed existing output folder. ") else: print("No existing output folder found. Nothing to remove.") else: print("Cleanup skipped. Set RUN_CLEANUP = True to remove the existing output folder.")

Subject: sub0045
Output folder: /Users/yc4421/Library/CloudStorage/Box-Box/GTA/ReCoDE_BrainMesh/ReCoDE-brain-mesh-creation/data/subjects/sub0045/output
No existing output folder found. Nothing to remove.

1. Create the initial mesh¶

The labelled geometry is converted into a mesh using a direct voxel-to-element approach. In simple terms, the labelled image is converted into a structured finite element mesh, where each region in the image becomes a corresponding part in the LS-DYNA keyword file.

The target mesh_size controls the approximate element size. A smaller mesh size gives a more detailed model, but the files become larger and simulations take longer. For the example subject, mesh_size = 2 is a practical default.

This step also performs voxel clean-up and creates additional anatomical structures used in the final head model, including the falx, tentorium, pia, and dura. These structures are thin membranes and are difficult to capture directly from MRI, so the pipeline reconstructs them algorithmically from surrounding anatomical landmarks.

In [8]:

# --------------------------------------------------
# Input and output files
# --------------------------------------------------
subject_path = SUBJECT_DIR
pre_model_path = PRE_MODEL_FILE
output_path = OUTPUT_DIR

print("=== Mesh creation ===")
print(f"Subject: {subject_name}")
print(f"Input geometry: {pre_model_path.name}")
print(f"Output folder: {output_path}")
print(f"Mesh size: {mesh_size} mm")
print()

# --------------------------------------------------
# Check input geometry
# --------------------------------------------------
if not subject_path.is_dir():
    raise FileNotFoundError(f"Cannot find subject folder: {subject_path}")

if not pre_model_path.is_file():
    raise FileNotFoundError(
        f"Cannot find pre_model.nii.gz at: {pre_model_path}, Please run 02_ImageProcess.ipynb first to generate the pre-model geometry."
    )

# The mesh creation code writes output files relative to the subject folder.
os.chdir(subject_path)
print(f"Working directory: {Path.cwd()}")
print()

# --------------------------------------------------
# [1/8] Initialise model creation
# --------------------------------------------------
print("[1/8] Initialising model creation")
print("      Input: pre_model.nii.gz")
print(f"      Output folder: {output_name}")

# create_model reads the labelled NIfTI geometry and prepares the output folder.
# DEPENDENCIES_DIR provides supporting files needed by the mesh creation code.
my_model = bmctk.create_model("pre_model.nii.gz", output_name, str(DEPENDENCIES_DIR))
print("      Status: model initialised")
print()

# --------------------------------------------------
# [2/8] Create coordinates
# --------------------------------------------------
print("[2/8] Creating model coordinates")

# This creates the coordinate system used to map image voxels into mesh nodes.
my_model.create_coordinates()
print("      Status: coordinates created")
print()

# --------------------------------------------------
# [3/8] Resample geometry
# --------------------------------------------------
print("[3/8] Resampling geometry")
print(f"      Target mesh size: {mesh_size} mm")

# Resampling controls the final element size. Smaller values retain more
# anatomical detail but increase the number of elements.
my_model.resample(mesh_size)
print("      Status: geometry resampled")
print()

# --------------------------------------------------
# [4/8] Clean tissue labels and repair local contacts
# --------------------------------------------------
print("[4/8] Cleaning tissue labels and local contacts")
print(f"      Voxel correction iterations: {iterations}")
print(f"      Neighbourhood blocks: {blocks}")
print(f"      Relabelling threshold: {threshold}")

# gm_check repairs grey-matter labels that are inconsistent or incomplete.
my_model.gm_check()

# check_contact(part_a, replacement_part, part_b) is used here to repair local
# interfaces where two tissues touch in a way that would be problematic for the
# mesh. The comments below describe the intended anatomical correction.

# Prevent skull from directly contacting cortical or cerebellar grey matter by
# restoring CSF between the brain surface and skull.
my_model.check_contact(257, 256, 42)  # Skull / CSF / Right cerebral cortex.
my_model.check_contact(257, 256, 3)   # Skull / CSF / Left cerebral cortex.
my_model.check_contact(257, 256, 47)  # Skull / CSF / Right cerebellum cortex.
my_model.check_contact(257, 256, 8)   # Skull / CSF / Left cerebellum cortex.

# Repair places where the CSF boundary meets background/void rather than skull.
my_model.check_contact(0, 257, 256)   # Void / skull / CSF.

# Prevent white matter from directly contacting CSF by restoring the grey matter
# layer where needed.
my_model.check_contact(41, 42, 256)   # Right white matter / right cortex / CSF.
my_model.check_contact(2, 3, 256)     # Left white matter / left cortex / CSF.
my_model.check_contact(46, 47, 256)   # Right cerebellar white matter / cortex / CSF.
my_model.check_contact(7, 8, 256)     # Left cerebellar white matter / cortex / CSF.

# Apply iterative voxel corrections to smooth out isolated or inconsistent labels.
my_model.voxel_corrections(iterations, threshold, blocks)
print("      Status: tissue labels corrected")
print()

# --------------------------------------------------
# [5/8] Create falx
# --------------------------------------------------
print("[5/8] Creating falx")

# The falx is a dural fold between the left and right cerebral hemispheres.
LR_center_coords = my_model.create_falx()
LR_center = LR_center_coords[0]
print(f"      Left-right centre coordinate: {LR_center}")
print("      Status: falx created")
print()

# --------------------------------------------------
# [6/8] Create tentorium
# --------------------------------------------------
print("[6/8] Creating tentorium")

# The tentorium is a dural fold separating the cerebrum from the cerebellum.
my_model.create_tentorium()
print("      Status: tentorium created")
print()

# --------------------------------------------------
# [7/8] Create LS-DYNA keyword mesh
# --------------------------------------------------
print("[7/8] Creating LS-DYNA keyword mesh")

# This writes the hexahedral mesh and part information to mesh.k.
max_node = my_model.create_k_mesh()
print(f"      Largest node number assigned: {max_node}")
print("      Status: mesh.k created")
print()

# --------------------------------------------------
# [8/8] Create node and element set list
# --------------------------------------------------
print("[8/8] Creating set list")

# The set list stores useful node and element groups for later visualisation,
# boundary conditions, contact definitions, or simulation setup.
my_model.set_list()
print("      Status: set list created")
print()

print("Mesh creation complete.")

# -------------------------------------------------- # Input and output files # -------------------------------------------------- subject_path = SUBJECT_DIR pre_model_path = PRE_MODEL_FILE output_path = OUTPUT_DIR print("=== Mesh creation ===") print(f"Subject: {subject_name}") print(f"Input geometry: {pre_model_path.name}") print(f"Output folder: {output_path}") print(f"Mesh size: {mesh_size} mm") print() # -------------------------------------------------- # Check input geometry # -------------------------------------------------- if not subject_path.is_dir(): raise FileNotFoundError(f"Cannot find subject folder: {subject_path}") if not pre_model_path.is_file(): raise FileNotFoundError( f"Cannot find pre_model.nii.gz at: {pre_model_path}, Please run 02_ImageProcess.ipynb first to generate the pre-model geometry." ) # The mesh creation code writes output files relative to the subject folder. os.chdir(subject_path) print(f"Working directory: {Path.cwd()}") print() # -------------------------------------------------- # [1/8] Initialise model creation # -------------------------------------------------- print("[1/8] Initialising model creation") print(" Input: pre_model.nii.gz") print(f" Output folder: {output_name}") # create_model reads the labelled NIfTI geometry and prepares the output folder. # DEPENDENCIES_DIR provides supporting files needed by the mesh creation code. my_model = bmctk.create_model("pre_model.nii.gz", output_name, str(DEPENDENCIES_DIR)) print(" Status: model initialised") print() # -------------------------------------------------- # [2/8] Create coordinates # -------------------------------------------------- print("[2/8] Creating model coordinates") # This creates the coordinate system used to map image voxels into mesh nodes. my_model.create_coordinates() print(" Status: coordinates created") print() # -------------------------------------------------- # [3/8] Resample geometry # -------------------------------------------------- print("[3/8] Resampling geometry") print(f" Target mesh size: {mesh_size} mm") # Resampling controls the final element size. Smaller values retain more # anatomical detail but increase the number of elements. my_model.resample(mesh_size) print(" Status: geometry resampled") print() # -------------------------------------------------- # [4/8] Clean tissue labels and repair local contacts # -------------------------------------------------- print("[4/8] Cleaning tissue labels and local contacts") print(f" Voxel correction iterations: {iterations}") print(f" Neighbourhood blocks: {blocks}") print(f" Relabelling threshold: {threshold}") # gm_check repairs grey-matter labels that are inconsistent or incomplete. my_model.gm_check() # check_contact(part_a, replacement_part, part_b) is used here to repair local # interfaces where two tissues touch in a way that would be problematic for the # mesh. The comments below describe the intended anatomical correction. # Prevent skull from directly contacting cortical or cerebellar grey matter by # restoring CSF between the brain surface and skull. my_model.check_contact(257, 256, 42) # Skull / CSF / Right cerebral cortex. my_model.check_contact(257, 256, 3) # Skull / CSF / Left cerebral cortex. my_model.check_contact(257, 256, 47) # Skull / CSF / Right cerebellum cortex. my_model.check_contact(257, 256, 8) # Skull / CSF / Left cerebellum cortex. # Repair places where the CSF boundary meets background/void rather than skull. my_model.check_contact(0, 257, 256) # Void / skull / CSF. # Prevent white matter from directly contacting CSF by restoring the grey matter # layer where needed. my_model.check_contact(41, 42, 256) # Right white matter / right cortex / CSF. my_model.check_contact(2, 3, 256) # Left white matter / left cortex / CSF. my_model.check_contact(46, 47, 256) # Right cerebellar white matter / cortex / CSF. my_model.check_contact(7, 8, 256) # Left cerebellar white matter / cortex / CSF. # Apply iterative voxel corrections to smooth out isolated or inconsistent labels. my_model.voxel_corrections(iterations, threshold, blocks) print(" Status: tissue labels corrected") print() # -------------------------------------------------- # [5/8] Create falx # -------------------------------------------------- print("[5/8] Creating falx") # The falx is a dural fold between the left and right cerebral hemispheres. LR_center_coords = my_model.create_falx() LR_center = LR_center_coords[0] print(f" Left-right centre coordinate: {LR_center}") print(" Status: falx created") print() # -------------------------------------------------- # [6/8] Create tentorium # -------------------------------------------------- print("[6/8] Creating tentorium") # The tentorium is a dural fold separating the cerebrum from the cerebellum. my_model.create_tentorium() print(" Status: tentorium created") print() # -------------------------------------------------- # [7/8] Create LS-DYNA keyword mesh # -------------------------------------------------- print("[7/8] Creating LS-DYNA keyword mesh") # This writes the hexahedral mesh and part information to mesh.k. max_node = my_model.create_k_mesh() print(f" Largest node number assigned: {max_node}") print(" Status: mesh.k created") print() # -------------------------------------------------- # [8/8] Create node and element set list # -------------------------------------------------- print("[8/8] Creating set list") # The set list stores useful node and element groups for later visualisation, # boundary conditions, contact definitions, or simulation setup. my_model.set_list() print(" Status: set list created") print() print("Mesh creation complete.")

=== Mesh creation ===
Subject: sub0045
Input geometry: pre_model.nii.gz
Output folder: /Users/yc4421/Library/CloudStorage/Box-Box/GTA/ReCoDE_BrainMesh/ReCoDE-brain-mesh-creation/data/subjects/sub0045/output
Mesh size: 2 mm

Working directory: /Users/yc4421/Library/CloudStorage/Box-Box/GTA/ReCoDE_BrainMesh/ReCoDE-brain-mesh-creation/data/subjects/sub0045

[1/8] Initialising model creation
      Input: pre_model.nii.gz
      Output folder: output
Reading brain data... 
      Status: model initialised

[2/8] Creating model coordinates
      Status: coordinates created

[3/8] Resampling geometry
      Target mesh size: 2 mm
      Status: geometry resampled

[4/8] Cleaning tissue labels and local contacts
      Voxel correction iterations: 3
      Neighbourhood blocks: 1
      Relabelling threshold: 0.9
Assigning missing brain matter to grey matter and brain stem... 
Number of missing voxels: 1033
Number of missing voxels: 887
Number of missing voxels: 847
Number of missing voxels: 833
Number of missing voxels: 827
Number of missing voxels: 822
Number of missing voxels: 822
Assigning remaining voxels... 
Number of missing voxels: 0
Voxel correction: 
Iteration 1
Values changed: 112

Voxel correction: 
Iteration 2
Values changed: 2

Voxel correction: 
Iteration 3
Values changed: 0

      Status: tissue labels corrected

[5/8] Creating falx
Falx location voxel: 45.51285702016037
Falx location scanner: 0.025714040320735876
      Left-right centre coordinate: 0.025714040320735876
      Status: falx created

[6/8] Creating tentorium
      Status: tentorium created

[7/8] Creating LS-DYNA keyword mesh
Reading nodal coordinates... 
Assigning element nodes... 
Writting K file...
      Largest node number assigned: 2964944
      Status: mesh.k created

[8/8] Creating set list
      Status: set list created

Mesh creation complete.

2. Smooth the mesh¶

The direct voxel-to-mesh conversion produces a blocky surface because the input geometry is voxel-based. This is expected, but sharp stair-step edges can create unrealistic local strain concentrations during simulation.

The smoothing step reduces these stair-step artefacts while trying to preserve the overall head and brain geometry. The result is saved as mesh_smoothed.k.

In [10]:

%%bash
set -euo pipefail

# --------------------------------------------------
# Input and output files
# --------------------------------------------------
subject_path="$SUBJECT_DIR"
output_path="$OUTPUT_DIR"
input_mesh="$OUTPUT_DIR/mesh.k"
output_mesh="$OUTPUT_DIR/mesh_smoothed.k"
smoother_path="$SMOOTHER_PATH"
mesh_size="$MESH_SIZE"

# --------------------------------------------------
# Smoothing parameters
# --------------------------------------------------
# These values are passed directly to the mesh_smoother executable.
# They control the smoothing strength, internal iteration/search limits, and
# the mesh-size setting used when checking the smoothed geometry. These defaults
# are the standard values used in this pipeline.
smoothing_passes=8
relaxation_factor=0.2
surface_iteration_limit=200
internal_search_limit=2000

echo "=== Mesh smoothing ==="
echo "Subject: $SUBJECT_NAME"
echo "Input mesh: ${input_mesh##*/}"
echo "Output mesh: ${output_mesh##*/}"
echo "Smoother: ${smoother_path##*/}"
echo

# --------------------------------------------------
# [1/3] Check required files
# --------------------------------------------------
echo "[1/3] Checking input files"

if [ ! -x "$smoother_path" ]; then
    echo "      Error: cannot find smoothing executable or it is not executable: $smoother_path"
    exit 1
fi

if [ ! -f "$input_mesh" ]; then
    echo "      Error: cannot find input mesh: $input_mesh"
    exit 1
fi

echo "      Status: required files found"
echo

# --------------------------------------------------
# [2/3] Move into mesh output folder
# --------------------------------------------------
echo "[2/3] Moving into mesh output folder"
echo "      Folder: $output_path"

cd "$output_path"
echo "      Working directory: $(pwd)"
echo

# --------------------------------------------------
# [3/3] Run mesh smoother
# --------------------------------------------------
echo "[3/3] Running mesh smoother"
echo "      Smoothing passes: $smoothing_passes"
echo "      Relaxation factor: $relaxation_factor"
echo "      Surface iteration limit: $surface_iteration_limit"
echo "      Internal search limit: $internal_search_limit"
echo "      Mesh size: $mesh_size mm"

"$smoother_path" mesh.k mesh_smoothed.k "$smoothing_passes" "$relaxation_factor" "$surface_iteration_limit" "$internal_search_limit" "$mesh_size"

[ -f "$output_mesh" ] && echo "      Status: smoothed mesh created"
echo

echo "Mesh smoothing complete."

%%bash set -euo pipefail # -------------------------------------------------- # Input and output files # -------------------------------------------------- subject_path="$SUBJECT_DIR" output_path="$OUTPUT_DIR" input_mesh="$OUTPUT_DIR/mesh.k" output_mesh="$OUTPUT_DIR/mesh_smoothed.k" smoother_path="$SMOOTHER_PATH" mesh_size="$MESH_SIZE" # -------------------------------------------------- # Smoothing parameters # -------------------------------------------------- # These values are passed directly to the mesh_smoother executable. # They control the smoothing strength, internal iteration/search limits, and # the mesh-size setting used when checking the smoothed geometry. These defaults # are the standard values used in this pipeline. smoothing_passes=8 relaxation_factor=0.2 surface_iteration_limit=200 internal_search_limit=2000 echo "=== Mesh smoothing ===" echo "Subject: $SUBJECT_NAME" echo "Input mesh: ${input_mesh##*/}" echo "Output mesh: ${output_mesh##*/}" echo "Smoother: ${smoother_path##*/}" echo # -------------------------------------------------- # [1/3] Check required files # -------------------------------------------------- echo "[1/3] Checking input files" if [ ! -x "$smoother_path" ]; then echo " Error: cannot find smoothing executable or it is not executable: $smoother_path" exit 1 fi if [ ! -f "$input_mesh" ]; then echo " Error: cannot find input mesh: $input_mesh" exit 1 fi echo " Status: required files found" echo # -------------------------------------------------- # [2/3] Move into mesh output folder # -------------------------------------------------- echo "[2/3] Moving into mesh output folder" echo " Folder: $output_path" cd "$output_path" echo " Working directory: $(pwd)" echo # -------------------------------------------------- # [3/3] Run mesh smoother # -------------------------------------------------- echo "[3/3] Running mesh smoother" echo " Smoothing passes: $smoothing_passes" echo " Relaxation factor: $relaxation_factor" echo " Surface iteration limit: $surface_iteration_limit" echo " Internal search limit: $internal_search_limit" echo " Mesh size: $mesh_size mm" "$smoother_path" mesh.k mesh_smoothed.k "$smoothing_passes" "$relaxation_factor" "$surface_iteration_limit" "$internal_search_limit" "$mesh_size" [ -f "$output_mesh" ] && echo " Status: smoothed mesh created" echo echo "Mesh smoothing complete."

=== Mesh smoothing ===
Subject: sub0045
Input mesh: mesh.k
Output mesh: mesh_smoothed.k
Smoother: mesh_smoother_macos_arm64

[1/3] Checking input files
      Status: required files found

[2/3] Moving into mesh output folder
      Folder: /Users/yc4421/Library/CloudStorage/Box-Box/GTA/ReCoDE_BrainMesh/ReCoDE-brain-mesh-creation/data/subjects/sub0045/output
      Working directory: /Users/yc4421/Library/CloudStorage/Box-Box/GTA/ReCoDE_BrainMesh/ReCoDE-brain-mesh-creation/data/subjects/sub0045/output

[3/3] Running mesh smoother
      Smoothing passes: 8
      Relaxation factor: 0.2
      Surface iteration limit: 200
      Internal search limit: 2000
      Mesh size: 2 mm
Dilatational wave speed [m/s] in: 
skin * cortical bone * diploe * CSF * brain * Falx * Tentorium * pia matter *
206.475 * 2873.84 * 1755.82 * 154.94 * 220.481 * 188.383 * 0 * 
Reading the input file ...

Number of solid elements: 370618
Number of shell elements: 87485
Number of nodes: 389919
Largest node number: 2964944

Reading the input file finished.

Node classification and neighbour listing ...

End of classification and neighbour listing.

Applying the smoothing filter ...

1 iterations done.
Minimum Jacobian (at element centre): 0.5904
Minimum characteristic length (2D - 3D): 1.48568 - 1.3608
Minimum time step: 0.000521809
2 iterations done.
Minimum Jacobian (at element centre): 0.353085
Minimum characteristic length (2D - 3D): 1.17299 - 0.991604
Minimum time step: 0.000413412
3 iterations done.
Minimum Jacobian (at element centre): 0.213474
Minimum characteristic length (2D - 3D): 0.958081 - 0.759081
Minimum time step: 0.000335346
4 iterations done.
Minimum Jacobian (at element centre): 0.130264
Minimum characteristic length (2D - 3D): 0.772318 - 0.602698
Minimum time step: 0.000273017
5 iterations done.
Minimum Jacobian (at element centre): 0.0801225
Minimum characteristic length (2D - 3D): 0.627494 - 0.491934
Minimum time step: 0.000226882
6 iterations done.
Minimum Jacobian (at element centre): 0.0496238
Minimum characteristic length (2D - 3D): 0.512738 - 0.410138
Minimum time step: 0.000182435
7 iterations done.
Minimum Jacobian (at element centre): 0.0309248
Minimum characteristic length (2D - 3D): 0.412415 - 0.347653
Minimum time step: 0.000177768
8 iterations done.
Minimum Jacobian (at element centre): 0.0193809
Minimum characteristic length (2D - 3D): 0.329355 - 0.298586
Minimum time step: 0.000177768

Applying the smoothing filter finished.

Writing the k file ...

DONE!
      Status: smoothed mesh created

Mesh smoothing complete.

3. Repair brain-skull contact¶

After smoothing, some brain surface nodes may become too close to the skull or may share contact with skull elements. This would be anatomically unrealistic because the brain should be separated from the skull by CSF and membranes.

This step checks the smoothed mesh and reclassifies problematic skull-contact elements as CSF where needed. The revised mesh is saved as mesh_smoothed_revised.k.

In [11]:

# --------------------------------------------------
# Input and output files
# --------------------------------------------------
input_k_file = OUTPUT_DIR / "mesh_smoothed.k"
output_k_file = OUTPUT_DIR / "mesh_smoothed_revised.k"

print("=== Brain-skull contact repair ===")
print(f"Subject: {subject_name}")
print(f"Input mesh: {input_k_file.name}")
print(f"Output mesh: {output_k_file.name}")
print()

# --------------------------------------------------
# [1/2] Check input mesh
# --------------------------------------------------
print("[1/2] Checking input mesh")

if not input_k_file.is_file():
    raise FileNotFoundError(f"Cannot find input mesh file: {input_k_file}")

print("      Status: input mesh found")
print()

# --------------------------------------------------
# [2/2] Create CSF buffer at the brain-skull interface
# --------------------------------------------------
print("[2/2] Creating CSF buffer")

# create_csf_buffer repairs non-physical brain-skull contact by restoring a CSF
# layer at the interface where required.
mesh_utils.create_csf_buffer(str(input_k_file), str(output_k_file))

if output_k_file.is_file():
    print("      Status: revised mesh created")
else:
    print("      Warning: command finished, but the revised mesh file was not found")
print()

print("Brain-skull contact repair complete.")

# -------------------------------------------------- # Input and output files # -------------------------------------------------- input_k_file = OUTPUT_DIR / "mesh_smoothed.k" output_k_file = OUTPUT_DIR / "mesh_smoothed_revised.k" print("=== Brain-skull contact repair ===") print(f"Subject: {subject_name}") print(f"Input mesh: {input_k_file.name}") print(f"Output mesh: {output_k_file.name}") print() # -------------------------------------------------- # [1/2] Check input mesh # -------------------------------------------------- print("[1/2] Checking input mesh") if not input_k_file.is_file(): raise FileNotFoundError(f"Cannot find input mesh file: {input_k_file}") print(" Status: input mesh found") print() # -------------------------------------------------- # [2/2] Create CSF buffer at the brain-skull interface # -------------------------------------------------- print("[2/2] Creating CSF buffer") # create_csf_buffer repairs non-physical brain-skull contact by restoring a CSF # layer at the interface where required. mesh_utils.create_csf_buffer(str(input_k_file), str(output_k_file)) if output_k_file.is_file(): print(" Status: revised mesh created") else: print(" Warning: command finished, but the revised mesh file was not found") print() print("Brain-skull contact repair complete.")

=== Brain-skull contact repair ===
Subject: sub0045
Input mesh: mesh_smoothed.k
Output mesh: mesh_smoothed_revised.k

[1/2] Checking input mesh
      Status: input mesh found

[2/2] Creating CSF buffer
Processing mesh: /Users/yc4421/Library/CloudStorage/Box-Box/GTA/ReCoDE_BrainMesh/ReCoDE-brain-mesh-creation/data/subjects/sub0045/output/mesh_smoothed.k
   -> Identifying Brain Nodes...
   -> Checking for Brain-Skull/Skin contact...
Writing result to: /Users/yc4421/Library/CloudStorage/Box-Box/GTA/ReCoDE_BrainMesh/ReCoDE-brain-mesh-creation/data/subjects/sub0045/output/mesh_smoothed_revised.k
Success. Total elements converted to CSF buffer: 25595
      Status: revised mesh created

Brain-skull contact repair complete.

Outputs¶

After running this notebook, the main mesh files should be:

data/subjects/<subject_name>/output/mesh.k

Raw hexahedral mesh created from the labelled geometry.

data/subjects/<subject_name>/output/mesh_smoothed.k

Mesh after surface smoothing.

data/subjects/<subject_name>/output/mesh_smoothed_revised.k

Final mesh after brain-skull contact repair.

This final revised mesh is the recommended file to inspect visually and use in later simulation setup steps.

References¶

• Betsholtz C, Engelhardt B, Koh GY, et al. (2024) Advances and controversies in meningeal biology. Nat Neurosci 27:2056–2072.
• Chen Y, Ostoja-Starzewski M (2010) MRI-based finite element modeling of head trauma: spherically focusing shear waves. Acta Mech 213:155–167.
• Darvishi V, Del Giovane M, David MCB, et al. (2025) Towards a deeper understanding of axonal injury mechanisms using a novel biomechanical model of normal pressure hydrocephalus. Journal details to be checked before publication.
• Giudice JS, Zeng W, Wu T, et al. (2019) An analytical review of the numerical methods used for finite element modeling of traumatic brain injury. Ann Biomed Eng 47:1855–1872.
• Hernandez F, Giordano C, Goubran M, et al. (2019) Lateral impacts correlate with falx cerebri displacement and corpus callosum trauma in sports-related concussions. Biomech Model Mechanobiol 18:631–649.
• Tadepalli SC, Erdemir A, Cavanagh PR (2011) Comparison of hexahedral and tetrahedral elements in finite element analysis of the foot and footwear. J Biomech 44:2337–2343.

In [ ]: