Tutorial 28 – Mesh Generation Based on Computer Tomography Data

Requirements:
• Partial Differential Equation Toolbox
• MATLAB implementation of the marching cubes algorithm, source
• stlwrite function for exporting to STL format, source

Introduction

This tutorial shows how to generate a mesh based on data obtained from a CT scan. Computer tomography is a relatively easy and straightforward way of reconstructing complex geometries, such as porous media. A demonstration of how to perform a simulation of fluid flow through a porous medium can be found in tutorial 29.

Data acquisition

First, we will demonstrate how to import .bmp images into MATLAB. We have 41 images, each 41×41
pixels, representing the scanned cube-shaped region.

number_of_images = 41;
image_width = 41;
image_height = 41;
file_extention = '.bmp';

The images are saved in the folder ‘images’.

image_folder = 'images';

We will now import them as logical values into the 41x41x41 array CT, which we first preallocate.

CT = false(41, 41, 41); % preallocate CT
for image_no = 1:number_of_images % Iterate in images
image_path = [image_folder, filesep, num2str(image_no), file_extention]; % images/image_no.bmp
CT(:, :, image_no) = logical(imread(image_path)); % import image as binary
end

Image processing

First, in order for the marching cubes algorithm to work correctly, we need to make sure the model is “watertight”, i.e. the “sides” of CT are set to false.

CT([1, end], :, : ) = false;
CT(:, [1, end], : ) = false;
CT(:, :, [1, end]) = false;

Next, we generate the grid at which the individual entries in CT were sampled.

x_spacing = linspace(.5, image_width - .5, 41);
y_spacing = linspace(.5, image_height - .5, 41);
z_spacing = linspace(.5, number_of_images - .5, 41);
[X, Y, Z] = meshgrid(x_spacing, y_spacing, z_spacing);

We now convert the geometry stored as voxels into the stl format using the marching cubes algorithm and then save it. For more information about the marching cubes algorithm the reader can refer to the paper by Lorensen and Cline it was originally proposed in .

[ matlab firstline=”17″][faces, vertices] = MarchingCubes(X, Y, Z, CT, .5);
model_filename = ‘porous_model.stl’;
stlwrite(model_filename, faces, vertices); [/matlab]

Mesh generation

We will use the PDE Toolbox for meshing our geometry. First, we initialize the model.

model = createpde();

We import the stl into the model.

importGeometry(model, model_filename);

Finally, we generate a linear mesh.

mesh = generateMesh(model, 'GeometricOrder', 'linear');

Converting the mesh from the PDE toolbox format to the QuickerSim CFD Toolbox format

First, we extract the p and t arrays from the model.

p = model.Mesh.Nodes;
t = model.Mesh.Elements;

Next, we scale the mesh.

p = p*1e-4;

We add the physical ID of the domain.

t(end + 1, : ) = 1;

We can now display the mesh.

displayMesh(p, t); To generate the e array containing boundary information, we will use the triangulation class available in MATLAB. We create an array containing the IDs of vertices lying on the boundary faces.

triangulation_object = triangulation(t(1:end - 1, :)', p');
faces = triangulation_object.freeBoundary';

We now have the faces of our mesh. We still need to assign them physical IDs. Let us first find the x coordinates of the face centers of mass.

face_centers_X = (p(1, faces(1, :)) + p(1, faces(2, :)) + p(1, faces(3, :)))/3;

The direction of the flow is aligned with the x axis, so the inlet faces will have the smallest x coordinate, and the outlet faces will have the largest x coordinates. We use logical indexing to translate this principle into code.

x_min = min(face_centers_X);
x_max = max(face_centers_X);
inlet_faces = face_centers_X < x_min + 1e-10;
outlet_faces = face_centers_X > x_max - 1e-10;

Note that we allow for a tolerance of 1e-8 due to the numerical accuracy of the meshing tool. We now create the rows of e corresponding to the faces’ physical ID and orientation.

e = [faces; ones(3, size(faces, 2))];
e(6, : ) = 0;

Finally, we assign the ID 1 to the walls, 2 to the inlet, and 3 to the outlet. We initially declared that all faces have ID 1, so we only need to set the IDs of the inlet and outlet.

e(4, inlet_faces) = 2;
e(4, outlet_faces) = 3;

We have now successfully imported the CT scan images, converted them into a geometric model, and meshed the model. We can now view it.

viewCase References

1. William E. Lorensen and Harvey E. Cline. 1987. Marching cubes: A high resolution 3D surface construction algorithm. SIGGRAPH Comput. Graph. 21, 4 (August 1987), 163-169. DOI: https://doi.org/10.1145/37402.37422

Author: Jakub Gałecki