2.1.

Overall Workflow

The overall workflow of this add-on is illustrated in Fig. 2. BlenderPhotonics contains three key steps. The first step is the construction of a surface mesh that is made of one or multiple objects. Users can create such surface models using the built-in shape objects provided by Blender or importing from user-defined data files as input. The second step is volume mesh generation and region labeling. The software exports the vertices and faces of the surface mesh generated from the first step and calls Iso2Mesh/TetGen to populate tetrahedra to fill the enclosed compartments of the surface mesh—each region is uniquely numbered. To permit mesh density and quality fine-tuning, a simple dialog is provided to allow for adjusting meshing parameters interactively, followed by regeneration of the volume mesh. The third step is mesh-based MC light simulation using the volume or regional meshes generated from the second step. The optical properties of each tissue region and the light source parameters are configured visually via the “custom property” in Blender, and an MMC compatible configuration data structure is generated and used as the input to launch the MMC simulation. The final light fluence distribution generated from such simulations, defined as floating-point fluence-rate values across the volume mesh, is then read in Blender and rendered for interactive visualization.

Fig. 2

Overall workflow diagram of BlenderPhotonics.

The modular software structure shown in Fig. 2 makes it effortless for future extensions. If desired, one can potentially replace the MMC simulation module with another mesh-based MC simulator, such as FullMonteCUDA.²⁴ Only lightweight code changes are required; these include adding additional “custom properties” in Blender, exporting these additional user inputs in the format acceptable by the target simulator, and importing the output data back to Blender for rendering.

2.2.

Installation

The installation of BlenderPhotonics requires setting up six software components: (1) Blender v2.8 or newer versions, (2) GNU Octave 4.2 or newer versions, (3) oct2py Python module (for interfacing Blender with GNU Octave), (4) jdata Python module (for reading/writing data exchange files),³⁷ (5) Iso2Mesh toolbox (for 3D mesh generation), and (6) MMCLAB toolbox (for photon MC simulation). Two optional components can be installed, including (1) ZMat toolbox (for reading/writing internally compressed mesh files) and (2) bjdata Python module [for reading binary JavaScript Object Notation (JSON)-based³⁷ data files]. All required software components are open source and are widely available. The step-by-step installation instructions can be found in Sec. 2.2 of the README file on the Github repository.³⁸

In our BlenderPhotonics package,³⁸ we also provide a fully automated installation shell-script for a Debian/Ubuntu-based Linux system. By simply replacing the apt command with yum or port, the attached installation script can be used to install BlenderPhotonics on Fedora Linux and Mac OS, respectively. For Windows users, although running the above shell script is possible if one has preinstalled Cygwin64, the apt command must be replaced because it is not supported on Windows. A Windows user may use winget install -e --id BlenderFoundation.Blender, winget install -e --id GNU.Octave, and winget install -e --id coti.mcxstudio in an administrator command window to install the needed software components, followed by installing the Python modules and Octave toolboxes manually.

Once all of the software components are successfully installed, one is able to start Blender, click on the menu “Edit/Preferences,” choose “Add-ons” in the Blender Preferences window, and search for “BlenderPhotonics”; the downloaded add-on should be listed (or use the “Install” button to browse a locally downloaded zip package). To enable BlenderPhotonics, one simply clicks on the check-box until a check-mark is shown. Although BlenderPhotonics has been tested on Blender 2.83, 2.92, and 3.0, the later versions are recommended because it supports an “exact” shape intersection solver starting from version 2.9. Aside from supporting GNU Octave, BlenderPhotonics also allows users to run Iso2Mesh and MMC via MATLAB (MathWorks, Natick, Massachusetts, USA) in the backend if one has installed the proprietary matlab.engine Python interface. One can choose between these two backends by a simple toggle button at the top of the BlenderPhotonics interface (see Fig. 3).

Fig. 3

BlenderPhotonics main interface and parameter dialogs (pointed by arrows) for each core feature.

2.3.

GUI Design

Ease-of-use is one of the key design requirements of BlenderPhotonics. To allow novice users to conveniently access all of the key functions of BlenderPhotonics, we created a simple function panel using Blender’s plug-in interface; see Fig. 3. In addition, it is desirable that a user is able to individually run any one of the functions. The independence of the task stages makes transferring or restarting a portion of the processing pipeline possible—e.g., using a portable computer to create the mesh model and then running GPU-based simulations on a dedicated server. The main interface allows for browsing volume data files, creating a volumetric mesh, importing regional mesh, and executing MMC simulation with only a single click.

2.4.

Data Inputs and Model Creation

BlenderPhotonics can create complex surface models (i.e., step 1 in Fig. 2) based on various forms of user inputs. Overall, there are three types of inputs that can be used for model creation. First, a user can directly add one or multiple built-in shape objects, such as cubes and spheres, using Blender and create a complex domain via Blender’s vastly rich built-in tools, such as Boolean modifier, bisect, knife tools, etc., and then create triangular/tetrahedral meshes via the “Blender2Mesh” module in BlenderPhotonics (see Fig. 3). Second, a user can import a triangular surface mesh pregenerated by other tools, handled by the “Surface2Mesh” module in BlenderPhotonics. Third, a 3D voxelated image, in the form of a Neuroimaging Informatics Technology Initiative (NIfTI),³⁹ JData³⁷/JNIfTI,⁴⁰ or a MATLAB/Octave .mat file, can be imported to create an image-derived surface mesh as supported in the “Volume2Mesh” module in BlenderPhotonics. Regardless of the input data type, the result of this step is a surface mesh defined by two data structures—vertices and faces. The vertices of the surface mesh are a floating-point array of dimensions $N_n\times 3$, with the three columns representing the $x/y/z$ coordinates, respectively, of a vertex in the surface mesh, where $N_n$ is the total number of nodes. The faces of the surface mesh are defined as a 2D integer array of size $N_e\times 3$, with each row containing three integer indices of the three nodes that form the triangular patch, where $N_e$ is the number of triangles. Only closed (i.e., watertight) surface mesh models are accepted because it is required for the subsequent tetrahedral mesh generation and MC simulation. A user should also pay attention to the complexity, measured by the number of nodes and triangles, of the surface model. Obviously, exceedingly dense surface meshes can lead to very dense volumetric mesh, resulting in excessively long meshing time and subsequent MC simulation run-times. A best practice is to create the surface mesh with the minimal number of nodes/triangles without losing the accuracy of the domain boundaries.

Blender is a vastly versatile environment for 3D shape modeling and domain creation, and how to use most of the functionalities of Blender for complex shape creation is beyond the scope of this tutorial. Readers should browse the large number of tutorials created by Blender’s developers and user community to learn how to effectively build complex 3D models using Blender. It is important to note that Blender separates built-in objects into (surface) mesh and nonmesh types, in which the former has well-defined node coordinates and face node indices and the latter may only have implicit surface definitions without explicitly defining surfaces or vertices. In many cases, nonmesh objects, such as a “metaball” object, can be converted to a mesh object by Blender during export. Other advanced built-in object types have not been tested.

A variety of surface mesh files can be imported to Blender. Blender supports most of the major 3D model files, such as OBJ geometry (.obj), STereoLithography (.stl), Filmbox (.fbx) formats, among others. Using BlenderPhotonics’s “Surface2Mesh” module and Iso2Mesh, we have also added support to import triangular surfaces stored in JSON/text-JMesh (.json or .jmsh), UBJSON/binary JMesh (.bmsh),⁴¹ object file (.off), and Infria MEDIT (.medit) formats. In addition, Blender, and subsequently BlenderPhotonics, can have extended file format support when a user installs appropriate add-ons to read/write a customized data format. For external model files, most models defined in these inputs contain the required vertex and face information. Similar to the model creation process, if there are objects of nonmesh types in the file, BlenderPhotonics automatically converts these objects to mesh types.

2.5.

Data Exchange between Blender and Octave via Portable JSON/JMesh Files

The data exchange between Blender and Octave is achieved via the combination of Blender’s Python API bpy, the Octave-Python interface oct2py, JSON-based data exchange Python module jdata, and the built-in JSON parser in Iso2Mesh. Although the default data exchange scheme in oct2py is achieved using MATLAB’s proprietary .mat format, we intentionally adopted the JSON-compatible JMesh format⁴¹ (.jmsh) because of the following advantages: (1) JMesh is human-readable whereas MATLAB.mat file is not, (2) JMesh is a plain JSON file that can be readily read/written in nearly all programming environments with lightweight parsers (many are built-in to the programming languages, such as Python, Perl, MATLAB, and JavaScript), (3) JMesh is directly editable and can be version-controlled whereas .mat files are binary only, and (4) last but not least, JMesh/JSON files can be readily used for web applications or hierarchical databases (such as NoSQL databases) whereas .mat files require specialized parsers that are not widely available. We want to highlight that JMesh is not a new format, but rather a set of JSON-compatible “name”/“value” pairs representing mesh-based data structures. Existing mesh data stored in other formats, such as Visualization Toolkit (VTK) format,⁴² can be potentially converted to JSON/JMesh and used in BlenderPhotonics via Octave-based parsers or external converters.

In Fig. 4, we use a simple tetrahedral mesh (of a unit cube) containing 8 nodes, 12 surface triangles, and 6 tetrahedral elements, shown in Fig. 4(a), as an example to illustrate the simplicity of using JSON/JMesh data annotations to store mesh data. Because the JMesh format is built upon the JData specification^37,43—a lightweight standard using JSON to annotate scientific data, JMesh mesh data constructs also support strongly-typed array [see Fig. 4(c)] and record-level binary data compression⁴³ [see Fig. 4(d)]. We strongly encourage developers of other biophotonics tools to consider adopting such a format in their applications to achieve easy interoperability between multiple tools and ease in future extensions.

Fig. 4

Sample JMesh representations of (a) a unit-cube (node numbering shown in circles). A JMesh file written with (b) plain JSON arrays, (c) JData-based annotated array, and (d) compressed binary array are supported and parsed by our Python and MATLAB toolboxes. Note that JMesh node indices (such as in “MeshTri3” and “MeshTet4” constructs) start from 1.

2.6.

Tetrahedral Mesh Generation from Blender Scenes

The surface mesh models created in Blender, consisting of a single or multiple objects, are used as the input for the next step—tetrahedral mesh generation. The mesh generation process consists of two parts: surface mesh preprocessing and tetrahedral mesh generation. The preprocessing of the surface mesh is achieved in Blender. First, Blender converts all nonmesh-type objects to mesh types to get well-defined points and faces. Other types of objects in the scene that cannot be converted into mesh objects, such as light sources and cameras, are excluded. All remaining objects are merged into one surface mesh object and checked for the presence of self-intersecting triangles. When self-intersection exists, Blender automatically inserts new points on the intersection line and split the faces, resulting in multiple subregions. At the end of the preprocessing step, a watertight nonself-intersecting triangular surface mesh is generated.

The preprocessed surface data are subsequently sent to Octave for mesh generation. Iso2Mesh is used to complete this step. First, the vertex coordinates and face information of the preprocessed model are saved as a .jmsh file and then loaded to Octave. In addition, several meshing parameters are prompt in a dialog (see Fig. 3) to the user to adjust the density of the output mesh, including (1) maximum tetrahedral element volume³⁶ and (2) percentage of surface mesh edges being kept after simplification, as defined in Iso2Mesh.³⁶ The latter parameter is passed to a mesh simplification algorithm implemented in the CGAL library^28,29 based on the Lindstrom–Turk algorithm.⁴⁴ The former parameter is passed to TetGen to set the maximum volume of the generated tetrahedral mesh. The two parameters act together to control the density of the generated mesh. In Octave, the inputs for Iso2Mesh are node coordinates, face data, and mesh generation parameters. In addition, Iso2Mesh calls TetGen (version 1.5) for mesh generation and region labeling. After tetrahedral mesh generation, Iso2Mesh returns three output data arrays, namely, “node,” “face,” and “elem” (see Fig. 5). They correspond to the vertices of the tetrahedral mesh, the faces triangle node indices, and the vertex indices of each tetrahedron, along with a label denoting the region to which it belongs, respectively. The tetrahedral mesh data are cached in a separate .jmsh file for subsequent calls. In addition, the exterior triangular surface for each tissue region/label is extracted and stored in a regional mesh .jmsh file to be loaded in Blender.

Fig. 5

Mesh data structure exchanged between Blender and Iso2Mesh.

The successful generation of a tetrahedral mesh readily enables a user to perform an array of advanced modeling tasks, such as performing finite-element (FE)- or boundary-element (BE)-based analyses to solve the forward or inverse problems in the domains such as mechanics, computational electromagnetics, computational optics, computational fluid dynamics, etc. The created volume and surface meshes can be exported from Blender in the format that is acceptable by specialized computational physics or multiphysics solvers, such as ANSYS, ABACUS, COMSOL, etc. Even if a user does not have the need to perform subsequent numerical analyses, Blender/BlenderPhotonics offers efficient visualization, manipulation, and advanced transformation to a surface or tetrahedral mesh, making rendering and processing complex 3D data easy.

2.7.

Tetrahedral Mesh Generation from 3D Volumetric Images

For researchers in the fields of medical imaging and biophotonics, creating simulations from 3D volumetric imaging data is a highly sought-after feature, but only limited tools are available for medical-image-based 3D mesh generation.²⁷ Iso2Mesh is one of only a few tools³⁶ that support this type of processing. In BlenderPhotonics, a dedicated submodule named “Volume2Mesh” is provided, as shown in Fig. 3. This module calls Iso2Mesh to automatically tessellate a 3D volume stored in the NIfTI,³⁹ JNIfTI,⁴⁰ or MATLAB.mat formats, and output a surface mesh model for subsequent processing. NIfTI is an open standard to store MRI/fMRI scans and is widely supported in the field of neuroimaging; JNIfTI is a standardized⁴⁰ JSON wrapper⁴³ to the NIfTI format to enable human-readability, easy extension, and data exchange between neuroimaging tools. By calling the streamlined mesh generation function in Iso2Mesh, the 3D volume data from a JNIfTI/NIfTI file can be preprocessed to create both tetrahedral volume mesh as well as surface mesh that is at the outer surface of the volume mesh. Multilabeled 3D volume can also be directly processed by Iso2Mesh to create the surface mesh model. Each labeled subregion is saved as a separate surface mesh record inside the exchange file. Unlike the previous two input types, the volume mesh and the regional mesh are generated simultaneously when a volume data file is imported into Blender. Nonetheless, a surface mesh extracted from the volume and regional mesh is displayed in Blender to allow for manual editing and further manipulation in case additional complex features, such as growing hairs as shown in the later section, are desired.

To use this feature, one simply clicks on the file browser button under the “Volume2Mesh” section of the BlenderPhotonics interface. This allows users to browse a text-JNIfTI (.jnii), binary-JNIfTI (.bnii), NIfTI (.nii/.nii.gz), or MATLAB.mat file stored in a local folder. In addition to loading a volume image from the user’s local disk, one can also type in a URL pointing to an online JNIfTI/NIfTI/.mat file. BlenderPhotonics automatically downloads the online data file and reads the content. After either a valid file path or URL is supplied, one can click on the “Convert 3D image file to mesh” button. A simple parameter dialog, shown in the middle of Fig. 3, pops up to allow users to set key meshing parameters, such as the upper-bound of tetrahedron volume ($V_{\max }$,^31,36 in cubic length unit), upper-bound of the Delaunay sphere radii of the surface triangles ($R_{\max }$,^28,36 in voxel unit), maximum allowed deviation (in voxel unit) from the voxelated boundary and mesh extract methods, including “cgalmesh,”²⁸ “cgasurf,”²⁹ and “simplify.”²⁷ Clicking on the “OK” button signals Iso2Mesh in Octave to start mesh generation. If successful, the regional mesh is loaded to Blender for inspection. If the mesh is not satisfactory, one can adjust the meshing setting and recreate the mesh.

2.8.

Mesh-Based Monte Carlo Photon Simulation Workflow

In this work, we are particularly interested in solving the RTE inside complex media using the MC method. Our previously developed MMC solver has already attracted a large and active user community consisting of students, researchers, and academics from biomedical optics and optical neuroimaging domains.²² However, the lack of an intuitive simulation domain preparation interface makes it difficult for less-experienced users to use. With the development of BlenderPhotonics, we specifically address this issue by interfacing the 3D mesh generation output from BlenderPhotonics with MMC and create a streamlined MC simulation environment in Blender. The MC simulation step requires three substeps: (1) domain preparation, (2) execution of photon simulation, and (3) output data visualization.

The domain preparation stage refers to the process of setting up the necessary optical simulation parameters, including optical properties for each tissue label, light source settings, and global simulation settings, such as total photon numbers, etc. After the completion of the mesh generation step above, Blender/BlenderPhotonics loads and visualizes regional meshes while assigning four custom properties for each region, i.e., ${\mu }_a$ (1/mm), ${\mu }_s$ (1/mm), $g$, and $n$ (Fig. 6). At the same time, a light source is added to the 3D domain. The light source has several user-customized properties, such as 3D position, orientation, source type, and total photon number to be simulated (Fig. 6). A user needs to set the optical parameters for each region of the mesh, place the light source in the desired position and adjust its orientation, and finally set the simulated photon number for the light source. It is important to note that the default unit of the domain is assumed to be mm in BlenderPhotonics. Once a user completes these setup tasks, he/she can then start the optical simulation using a dedicated button on the BlenderPhotonics panel.

Fig. 6

Optical parameters and light source configuration interface. The panels for setting optical parameters for region 1 (right) and region 2 (middle) are shown at the bottom; the panel for setting the source is shown at the middle-top.

When a photon simulation starts, BlenderPhotonics automatically reads the optical parameters of all-region models in Blender and the light source information and passes it into Octave. Octave automatically generates a configuration file for the light simulation based on the supplied data. The mesh data of the simulation domain is derived from the mesh data saved during the tetrahedral generation stage. Other critical simulation parameters are passed from the user’s settings during the preparation stage. In particular, the light source direction is stored in Blender as a quaternion. In Octave, the direction in quaternions is used to compute the standard vector form.⁴⁵ The total run-time of the MMC simulation depends on the model complexity, the total number of photons to be simulated, and the optical properties. A progress bar, alongside other simulation-related information, is printed in a command-line window when MMC is being executed.

After the MMC simulation is completed, the fluence map computed over the 3D mesh is saved into a JMesh file and loaded back to Blender for visualization and postprocessing. The fluence map is converted in the log-10 scale for better rendering of the simulation results. The log-scaled fluence intensity map is assigned as the “vertex weights” in Blender and rendered in pseudocolors using a built-in color map.

2.9.

BlenderPhotonics File Structure and Optimization

BlenderPhotonics interacts with three types of files—Python scripts, Octave scripts, and intermediate files. All Python scripts are placed in the add-on’s main directory. Octave scripts are stored in the “script” subfolder inside the add-on’s directory. All intermediate outputs and data exchange files (JMesh files) are stored inside a temporary directory—for Linux/MacOS, this folder is typically /tmp/iso2mesh-USER/blenderphotonics, and for Windows, this is typically C:\Users\USER\AppData\Local\Temp\iso2mesh-USER\blenderphotonics, where “USER” is replaced by the actual user account name.

In terms of run-times, BlenderPhotonics does not have a direct impact on either those of Iso2Mesh mesh generation or MMC light simulation because they were performed in Octave. However, when writing and reading large numbers of surface mesh files, there is a noticeable data transfer overhead. Therefore, when loading the mesh data to Blender for inspection or rendering MMC simulation results, we only extract the surfaces of the regional mesh to reduce such overhead.

3.1.

“SkinVessel” Benchmark—Creating Simulations from Blender Objects

In the first example, we show mesh generation and light simulation of multilayered skin tissues with an embedded blood vessel, adapted from the “SkinVessel” benchmark initially created by Dr. Steven Jacques for his MC software mcxyz.⁴⁶ The model consists of a multilayered slab structure, derived from a combination of a cube, dissecting planes, and a blood vessel created from a cylinder object in Blender; see Fig. 7. From bottom to top, the layers are named “Low-slab,” “Mid-slab,” and “High-slab,” respectively. The geometric and optical parameters of each region are reported in Table 1. The construction process of the model is shown in Figs. 7(a)–7(c).

Fig. 7

Intermediate steps of creating the SkinVessel benchmark: (a) after adding Blender objects, (b) after clicking “preview surface tesselation,” (c) output regional mesh, and (d) fluence rate cross-section generated by MMC photon simulation.

Table 1

Geometric and optical parameters of the “SkinVessel” benchmark.

To create this domain, we first use the menu “Add/Mesh” in the object-mode, to add a “Cube,” two “Plane” objects (plane-1 and plane-2), and a “Cylinder.” For each object, in the popup property dialog, we set their sizes (depth for the cylinder), offsets $x_0,y_0,z_0$, rotation angles ${\varphi }_x,{\varphi }_y,{\varphi }_z$ relative to each axis, and other properties according to Table 2. Here, we set all objects lengths in voxel unit to match the original benchmark made for mcxyz; the voxel size (in mm) is specified by the “unitinmm” property shown in Table 2. This results in a three-layered structure similar to Fig. 7(a). Once this model is created, one can click on the “Preview surface tesselation” button in BlenderPhotonics’s interface and inspect the tessellated surface mesh, as shown in Fig. 7(b). One should see no self-intersecting triangles, and each compartment of the surface must be watertight. After verification, select the menu “Edit/Undo” to return to the polyhedral surface model. In the next step, one clicks on the “Convert scene to tetra mesh” button on the BlenderPhotonics interface. In the shown dialog, we set the maximum tetrahedron volume to 30 (cubic length unit), optionally uncheck the “Convert to triangular mesh” while leaving other settings to default values. Clicking on “OK” allows Blender to call Iso2mesh and TetGen to create a tetrahedral mesh of 109,110 nodes and 623,608 tetrahedra. Once completed, the volume mesh is loaded into the domain for inspection. This usually takes about 5 to 10 s.

Table 2

Blender properties set for each object to build the “SkinVessel” simulation.

The next step is to configure an MMC simulation domain. To do so, one clicks on the “Load mesh and setup simulation” button in the “Multiphysics Simulations” section of BlenderPhotonics’s panel. This reloads the generated mesh in the previous step in the form of a regional mesh (instead of a volumetric mesh) and names each tissue region as “region_i” in the object list [Fig. 7(c)]. It also attaches default optical properties to each regional surface. In addition, a “light” object named “source” is also added to the scene positioned above the domain’s center. The custom property setting panels for each object are similar to those in Fig. 6. One should set the optical and source properties following Table 1. In this example, we use the settings for “source-1.” Once completed, the domain is ready for the next step—MMC simulation—which is achieved conveniently by a single click on the “Run MMC photon simulating” button. In the pop-up dialog, one can adjust the simulation setting, such as maximum time-gate, reflection settings, normalization, or choosing different GPU devices according to MMC’s command-line options. Clicking on the “OK” button starts the MMC simulation in the background. Once the simulation is completed, the computed light fluence map is loaded to the Blender rendering window and displays as “vertex weight” using a default color map [Fig. 7(d)].

3.2.

“Colin27” Benchmark—Creating Simulations from Surface Meshes

In the second example, we demonstrate the steps needed to create tetrahedral mesh models and run MMC simulations using precreated surface mesh models. The anatomical model was derived from a widely used human brain atlas, known as the “Colin27” atlas. A set of previously generated²² triangular surface meshes containing four closed tissue surfaces—scalp, CSF, gray matter (GM), and white matter (WM)—are saved in the form of a JSON/JMesh file and loaded to Blender using the “Import surface mesh” button of BlenderPhotonics; a cropped input surface mesh can be found in Fig. 8(a). By a single click on the “Convert scene to tetra mesh” button in BlenderPhotonics and setting 100 as the maximum element size, one can create a tetrahedral mesh with a cross-sectional plot shown in Fig. 8(b). Next, clicking on the “Load mesh and setup simulation” button, one can set the displayed regional meshes [Fig. 8(c)] using the optical parameters listed in Table 3, with the light source position, light source direction, and photon number being [75.76, 66.99, 168.21] mm, [0.1636, 0.4569, −0.873] mm, and ${10}^8$, respectively. Finally, a single click on the “Run MMC photon simulation” button starts the MMC simulation using the input atlas and loads the computed fluence map back to Blender once the simulation is completed [see Fig. 8(d)].

Fig. 8

Intermediate steps of creating the Colin27 benchmark: (a) cropped input surface mesh, (b) cropped tetrahedral mesh, (c) region mesh and source, and (d) cropped fluence overlaid on the regional mesh.

Table 3

Optical properties for each tissue type in the Colin27 benchmark.

3.3.

“Digimouse” Benchmark—Creating Simulations from Segmented Volumetric Images

In this example, we show the processing steps to convert a 3D volume into a mesh model and subsequently run an MC simulation. The benchmark is derived from a public dataset known as the “Digimouse” atlas, a segmented CT image of size $190\times 496\times 104\text{voxels}$ with an isotropic voxel size of 0.8 mm. The original CT images of the Digimouse atlas are converted to a JSON/JNIfTI file. To reproduce this example, one can directly type the file link⁴⁷ in the JNIfTI file field in BlenderPhotonics’s interface. One can also download this file and browse it on the local disk. Next, one clicks on the “Convert 3D image file to mesh” button and sets the maximum element volume to 100, and the volume image is loaded and processed in Iso2Mesh to create a tetrahedral mesh¹¹ made of 21 tissue types [Fig. 9(b)]. The regional mesh for each tissue type [Fig. 9(a)] is then loaded back to Blender, and the optical parameters of each region are assigned manually in the Blender property dialog, according to those listed in Table 4. The light source position, light source direction, and photon number are set to [40, 160, 80] mm, [0, 0, −1] mm, and ${10}^8$, respectively. One also sets the “unitinmm” property of the source to the input voxel size, which is 0.8 (mm) in this case. The tetrahedral mesh of the atlas consists of 72,815 vertices and 407,739 tetrahedra. The final MMC photon simulation results are shown in Figs. 9(c) and 9(d). To demonstrate the advanced rendering capability of Blender, we show two ways of creating cross-sectional images of the output. In Fig. 9(c), we first box-select a region of faces in the edit-mode and press “delete” on the keyboard. This creates a trimmed mesh showing the internal surfaces with the associated vertex-weight (fluence values). However, to create a flat cross-section, one can select “Mesh/Bisect” in the edit-mode and draw a straight line across the mesh. In the “Bisect” setting dialog, one can choose to remove either side of the bisecting plane and fill the cross-sectional plane with flat patches. The result is shown in Fig. 9(d). An even faster way to crop a mesh is to use the view-port cropping feature by pressing shortcut “Alt+B” and box-select the region to be kept.

Fig. 9

Intermediate steps of creating the Digimouse benchmark: (a) regional mesh after converting the volume to mesh, (b) output tetrahedral mesh (cropped), (c) fluence-rate map rendering over a trimmed mesh using face-deletion, and (d) that rendered over a sliced mesh using the “bisect” function.

Table 4

Optical parameters for each of the regions/labels of the “Digimouse” benchmark.48,49

3.4.

Advanced Modeling Example 1—Simulating Rough Surfaces Using BlenderPhotonics

BlenderPhotonics’ ability to bridge Blender’s superior model creation/transformation capability with state-of-the-art 3D mesh generation and quantitative photon simulation opens numerous possibilities to advance our understanding of complex tissue-photon iterations in a realistic setting, helping us design better instruments and experiments and accounting for scenarios that simple models could never reveal. We can not enumerate all possible advanced modeling functions provided by Blender; however, just to provide inspiration for biomedical optics researchers and computational scientists, we hand-picked two features that are relevant to tissue optics. In the first example, we use the rough-surface creation feature in Blender and investigate how changing the roughness of a skin-mimicking surface could lead to different light distributions and results. In the second example, we combine BlenderPhotonics with one of our most recent advances in MC simulation—the implicit MMC (iMMC) algorithm²⁶—to potentially enable the study of the impact of human hairs in fNIRS measurements. To the best of our knowledge, these types of studies have not been reported by other MC simulators.

In the rough-surface simulations, we adjusted the geometric parameters of the interface based the aforementioned “SkinVessel” benchmark to simulate the roughness of realistic skin surface, which was difficult to achieve in the past due to the lack of relevant features in the mesh generation tool. In this case, all parameters used in this benchmark are largely the same as we described in the previous subsection, except that we change the source position from those for “source-1” to those of “source-2” according to Table 2. We also shrink the domain by a factor of 10 to amplify the effects by changing the length unit, defined as a custom property named “unitinmm” attached to the source object, from 0.005 to 0.0005 mm. In addition, we modified the plane-1 object (see Table 1) to create a rough surface. This is achieved by first selecting the plane object in the edit mode, right-clicking on the object, choosing “Subdivide,” and setting the division number to 39 [Fig. 10(a)]. In the next step, one switches to the vertex-selection mode, selects the menu “Mesh/Transform/Randomize,” and then sets “Amount” to 1 and “Normal” to 1.0. This creates a maximum $\pm 1\times 0.0005\mathrm{mm}$ normal-direction-constrained random movement [Fig. 10(b)]. According to the literature,⁵⁰ the roughness of the human skin surface can be described using the parameter $R_a$, defined as the arithmetic mean deviation of the depth profile

Fig. 10

Comparisons between (b, e) low-roughness and (c, f) high-roughness surfaces using BlenderPhotonics: (a) initial surface model after subdivision, (b, c) created rough surfaces using randomized nodal offsets, (d) regional mesh, and (e, f) MMC simulated fluence rate cross-sectional plots from a tilted disk source.

The images shown in Figs. 10(e) and 10(f) demonstrate that considering tissue surface roughness can noticeably alter the incident light direction. This is because a rough surface acts as a diffuser and broadens the incident beam, significantly reducing the collimation of the incident beam. The capability of BlenderPhotonics to simulate such a realistic surface model can assist researchers in better optimizing their imaging system and obtaining more accurate quantification.

3.5.

Advanced Modeling Example 2—Simulating Realistic Human Hairs Using Blender Hair System

In addition to the roughness of realistic human skin, hair is also widely considered an important factor when performing optical measurements on the human body. In particular, the presence of human hair is considered one of the most challenging aspects of fNIRS measurements.⁵¹ Unfortunately, due to the high complexity of modeling realistic hair, previous publications on hair modeling were limited to simple models involving only dozens of hairs of uniform hair-root distribution and tilting angle. Here, we combine the advanced hair modeling system in Blender with the latest advances in the implicit MMC (iMMC) simulation algorithm and provide a viable path for researchers to rigorously simulate the effects of hair in realistic measurements. A comprehensive study characterizing the impact of hair of different densities, colors, lengths, etc., in the context of fNIRS will be reported in a separate publication. Here, we just want to demonstrate the rich and advanced hair model generation capability provided by Blender and hopefully motivate readers to explore other complex shape modeling features provided by this open-source platform.

In this example, we first show results from “growing” realistic hairs over a simple three-layered head model and then show hair models created using the aforementioned Colin27 atlas. First, we create a three-layered slab model similar to the steps used in Sec. 3.1. The slab has dimensions of $20\times 50\times 34.51{\mathrm{mm}}^3$, a 12.44-mm top layer simulating scalp/skull, a 2.07-mm middle layer simulating CSF, and a 20-mm bottom layer simulating the brain (GM/WM). To grow hairs on the top surface, we first select the top surface in the edit mode, right click and select “Subdivision,” with the number of cuts set to 1. Then we select the face-center-vertex inserted by the “Subdivision’’ step and open the “Object properties” panel on the right. Under the “Vertex group,” add a new group, and click the “Assign” and “Select” buttons to apply particle simulations over the selected surface later. Next, one switches to the “Particle properties” panel, clicks on the plus sign to create a new particle property, and then selects “Hair.” This setting immediately grows 1000 straight (along the normal direction) hairs with random roots on the entire surface of the slab. To constrain the hairs to only the top surface, go to the “Vertex group” subpanel, under the particle properties, and choose the vertex group created earlier. Once a user switches to the object model, the hairs should be displayed similar to Fig. 11(a). One can use the “Particle properties” panel to adjust various properties of the hair model, such as total count, length, tilting angle, randomness of the tilting angle, etc. For example, setting the $x/y/z$ velocity under the “Velocity” panel changes the direction of the hairs and randomness of the tilting angle; setting hairs to tilt 45 deg in the $+x$-axis with a randomness of 0.2 results in the hair model shown in Fig. 11(b). To create even more realistic hair shapes, one can check the “Hair Dynamics” checkbox and click the “Play” button at the bottom of the window. This applies forces to each segment of the hair shaft and creates curved hairs pulled by gravity. Dynamically curved hairs using Fig. 11(b) as the initial model results in Fig. 11(c). Similarly, this process can also be repeated in arbitrarily complex shapes, such as the Colin27 atlas. In Figs. 11(d) and 11(e), we show the outcome of growing 80,000 hairs ($165\text{hairs}/{\mathrm{cm}}^2$) with and without hair dynamics, respectively, on the Colin27 head surface; in Fig. 11(f), we increase the hair count to 400,000 ($825\text{hairs}/{\mathrm{cm}}^2$) with hair dynamics enabled to simulate extremely dense and realistic human hair. A special script is used to export all hair vertices and the head mesh into a file that is processed by the Iso2Mesh toolbox to create iMMC simulation models for subsequent light transport simulations.²⁶

Fig. 11

Demonstrations of different hair models created in Blender: (a) straight, (b) 45-deg tilted, and (c) gravity-curved hairs grown over a three-layered slab head model and (d) straight and (e, f) gravity-curved hairs simulated over the Colin27 atlas at two densities of (d and e) $165\text{hairs}/{\mathrm{cm}}^2$ and (f) $825\text{hairs}/{\mathrm{cm}}^2$.