This page illustrates a typical MIRAN workflow for going from a CT 3D image to a finite element model capturing the anatomy of the thorax. The goal is to illustrate MIRAN main functionalities and the typical steps involved in this process. The workflow demonstrated below is oriented towards building a thorax model for Electrical Impedance Tomography (EIT). The complete workflow in practice takes about 10 minutes. Python code snippets are included which demonstrate the process step-by-step, as well as a full Python script (including necessary the necessary Python modules) to download here.
CT DICOM Dataset
As an demonstration dataset we used the “ARTIFIX A Aorta w-c 1.5 B20f 60%” series from the ARTIFIX dataset.
The dataset contains CT images of the thorax, two slices from the 3D CT volume are illustrated below:
im = nes_dcm.LoadDicomVolume("ARTIFIX A Aorta w-c 1.5 B20f 60%", file_extension="*.dcm", study_number=0, series_number=0)
Surface Body Capture
The surface of the body is captured automatically by MIRAN using algorithms based on CT pixel value and connectivity information as illustrated in the images below. The extracted surface is “pixelated” meaning that the surface is not smooth as the original resolution of the CT image results in pixels that are visible on the surface. The extracted body surface also includes parts of the CT bed, which are later removed automatically.
Surface Body Smoothing
A data driven smoothing algorithm is used to remove the effect of pixels on the surface of the body. The strength of the smoothing effect can be controlled by the user.
body_surf_high_res = nes_msh.ExtractBodySurfaceFromCtImage(im, ct_value_threshold=-200, smoothing_steps=50, display_intermediate_steps=False)
Body Surface Decimation
The extracted body surface is generally composed of a large number of triangles, in the example above about 2.5 millions. This number of triangles is too large if the surface is to be fed into a volume meshing algorithm. A surface decimation function based on the ACVD algorithm has been incorporated in MIRAN, allowing effective reduction of the number of surface triangles while preserving the surface anatomical information.
In the above surfaces the number of triangles has been reduced from 2.5 millions to 19,000. The final number of triangles in the decimated mesh is controllable by the user.
body_surf_low_res = nes_msh.DecimateSurface(body_surf_high_res, num_output_points=9000)
CT Table Removal
CT tables contain metallic parts that are radio-opaque. Segmentation algorithms that are based on pixel intensity values and thresholds (like the algorithm used for extracting the body surface) are prone to segment also these parts of the table. The table is later removed automatically.
The images above show the result from using the automatic CT table removal algorithm. The image on the left shows the front of the body and the image on the right the back of the body.
body_surf_no_table = nes_msh.KeepLargestConnectedComponent(body_surf_low_res)
wr = vtk.vtkDataSetWriter()
wr.SetFileTypeToBinary()
wr.SetInputData(body_surf_no_table)
wr.SetFileName("BodySurfaceLowResNoTable.vtk")
wr.Write()
Segmentation of Lungs
Lungs appear in CT images as low intensity regions of the body (with a typical Hounsfield value of -500, compared to -50 for fat, and +40 for muscle). Lungs are well separated therefore in terms of CT image value from the rest of the organs and are easily segmented. MIRAN provides Level Set algorithms that implement the confidence connected criterion and which are based on a modified and enhanced version of the Sparse Level Set algorithms developed by Shawn Lankton.
Algorithm Initialization
The segmentation is initialized by placing a ball inside the organ of interest, this can be done graphically in the graphical version of MIRAN, or by specifying the ball coordinates and radius in the scripting version of MIRAN. Based on this initial ball the algorithm establishes the mean values and standard deviation of the pixels included in the ball and starts adding neighboring pixels that have similar pixel intensity statistics to the initial volume represented by the ball. The boundary of the organ is found when all the neighboring pixels do not match the intensity statistics of the initial volume, and therefore the growth of the segmented volume stops.
#initialization point in right lung
right_lung_point = (73,144,86)
#initialization sphere vtkFloarArray of (x,y,z,r) params
initialization_sphere_r = vtk.vtkFloatArray()
initialization_sphere_r.SetNumberOfComponents(4)
initialization_sphere_r.InsertNextTuple4(right_lung_point[0],right_lung_point[1],right_lung_point[2],15)
#segment right lung
right_lung_surface = nes_msh.SegmentAnOrgan(im, initialization_sphere_r, ls_iter=5000, ls_std_mult=1.5, ls_lambda=1.5, smoothing_iter=50, display_intermediate_steps=True)
#write to vtk format
wr = vtk.vtkDataSetWriter()
wr.SetFileTypeToBinary()
wr.SetInputData(right_lung_surface)
wr.SetFileName("SegmentedRightLung.vtk")
wr.Write()
#decimate right lung surface
right_lung_surface_low_res = nes_msh.DecimateSurface(right_lung_surface, num_output_points=4000)
#initialization point in left lung
left_lung_point = (219,135,86)
#initialization sphere vtkFloarArray of (x,y,z,r) params
initialization_sphere_l = vtk.vtkFloatArray()
initialization_sphere_l.SetNumberOfComponents(4)
initialization_sphere_l.InsertNextTuple4(left_lung_point[0],left_lung_point[1],left_lung_point[2],15) # 15 = radius in mm
#segment right lung
left_lung_surface = nes_msh.SegmentAnOrgan(im, initialization_sphere_l, ls_iter=5000, ls_std_mult=1.5, ls_lambda=1.5, smoothing_iter=50, display_intermediate_steps=True)
#write to vtk format
wr = vtk.vtkDataSetWriter()
wr.SetFileTypeToBinary()
wr.SetInputData(left_lung_surface)
wr.SetFileName("SegmentedLeftLung.vtk")
wr.Write()
#decimate left lung surface
left_lung_surface_low_res = nes_msh.DecimateSurface(left_lung_surface, num_output_points=4000)
Segmentation Results
The result of the segmentation, for the two lungs, is illustrated in the renderings below. Two initialization balls, one for each lung, have grown to capture the surface of the lungs. The algorithm includes smoothing functions and surface decimation functions as for the body surface capture illustrated above.
v1 = nes_viz.DisplayDataSetSurface(body_surf_no_table, color=(1,1,1), edge_color=(1,1,1), opacity=0.4, wireframe=True) v2 = nes_viz.DisplayDataSetSurface(right_lung_surface, color=nes_msh.lung_color, renderer=v1.ren) v3 = nes_viz.DisplayDataSetSurface(left_lung_surface, color=nes_msh.lung_color, renderer=v1.ren) nes_viz.RunRenderWindowInteractorOnDisplayedObject(v1)
Electrode Coordinate Generation
MIRAN offers functions for adding electrodes to the body surface, allowing simulating the application and sensing of electric potentials on the body, as occurs in Electrical Impedance Tomography and in other biomedical applications. A first set of these functions allows generating electrode coordinates, for example by intersecting the body with a plane and by generating points on the intersection line at equal distances around the body.
In the above illustration a semi-translucent intersection plane is rendered. Sixteen equispaced points are generated on the line resulting from the intersection of the body and the plane, these points are indicated as small green spheres.
#z level for each electrode plane, no. electrode per plane
planes = [60, 140]
num_elecs_per_plane = 16
num_planes = len(planes)
num_elecs = num_elecs_per_plane*num_planes
#storage for electrode coordinates
elec_coords = np.zeros((num_elecs_per_plane*num_planes,3), dtype=np.float32)
bounds = body_surf_no_table.GetBounds()
mid_x = (bounds[0]+bounds[1])/2
mid_y = (bounds[2]+bounds[3])/2
# for each z level generate equispaced electrodes on the specific z level
for i in range(num_planes):
# (0,0,1) is normal to plane used for intersection, here normal along z
intersect = nes_msh.FindEquispacedElectrodeCoordinatesAlongIntersectionPlane(body_surf_no_table, (mid_x, mid_y, planes[i]), (0,0,1))
# intersect is an object on which methods can be called
# finds length of polyline resulting from intersection of plane with body, and generate equipspaced electrode
length = intersect.GetIntersectionLength()
len_step = length / num_elecs_per_plane
# parametrized positions for electrodes we want to generate
np_lengths = np.arange(0,length,len_step)
list_of_lengths = list(np_lengths)
# use again the intersect object to generate (x,y,z) coordinates given the parametrized positions of the electrodes
elec_coords[i*num_elecs_per_plane:(i+1)*num_elecs_per_plane,:] = intersect.FindElecXYZCoordsFromParamterizedLength(list_of_lengths, display_intermediate_steps=True)
# given electrode coordinates we need to generate surface patches (round electrodes) on body surface and embed internal organ surfaces into a single vtkPolyData object for meshing
surf_to_be_mesehd = nes_msh.AddElectrodePatchesToBodySurfaceAndJoinResultWithInternalOrganSurfaces(elec_coords, body_surf_no_table, internal_surfaces_list=[right_lung_surface_low_res, left_lung_surface_low_res])
Volume Meshing and Generation of Electrodes on the Mesh Surface
Based on the extracted body surface, on the surfaces of the internal organs of interest (lungs in this example) and on the generated electrode coordinates MIRAN generates tetrahedral volume meshes that:
- Conform to the body surface.
- Conform to the internal organs (i.e. the mesh has internal sub-regions that match the internal organs geometry).
- Present surface patches that match the electrode positions and sizes.
- Have a specified mesh density in the different sub-regions.
- Have a specified (usually high) mesh density under the patches describing the electrodes.
MIRAN uses the TetGen mesh generation library and custom functions to achieve the above goals. Below is illustrated the surface of a produced volume mesh with two rings of 16 electrodes each. The mesh edges are rendered in white color to provide a visual representation of the mesh density.
tet = nes.vtkNESTetrahedralizePolyData()
tet.SetInputData(surf_to_be_mesehd)
# specify facet area constraints
facetConstraintIndex = vtk.vtkIntArray()
facetConstraintIndex.SetNumberOfComponents(1)
facetConstraintIndex.SetNumberOfTuples(num_elecs)
#constraint 0 applies to elec #1, constraint 1 to elec #2 etc
for i in range(num_elecs):
facetConstraintIndex.SetTuple1(i,i+1)
facetConstraintValue = vtk.vtkDoubleArray()
facetConstraintValue.SetNumberOfComponents(1)
facetConstraintValue.SetNumberOfValues(num_elecs)
# specifies for each constraint the maximum area
for i in range(num_elecs):
facetConstraintValue.SetTuple1(i,0.5)
tet.SetFacetConstraints(facetConstraintIndex, facetConstraintValue)
# point outside of lungs and inside body for defining the body region
body_point = (150,141,295)
# each region is specified with 5 parameters (x,y,z, subvolume_index, subvolume_volume_constraint)
regions = vtk.vtkDoubleArray()
regions.SetNumberOfComponents(6)
regions.SetNumberOfTuples(3)
regions.SetTuple6(0,right_lung_point[0],right_lung_point[1],right_lung_point[2],2,40,0) # first param is tuple index, last 0 argument is unused
regions.SetTuple6(1,left_lung_point[0],left_lung_point[1],left_lung_point[2],3,40,0)
regions.SetTuple6(2,body_point[0],body_point[1],body_point[2],1,120,0)
tet.SetVolumeRegionsInfo(regions)
# a = enforce volume constraints, A = apply tetrahedra attributes, p = mesh PLC, q = quality mesh, V = verbose , Q = quiet
tet.SetTetGenOptionsString("aApqQ")
tet.Update() # perform meshing
# extract volume and face elements
mesh = tet.GetOutputTetraMesh()
faces = tet.GetOutputTriFaces()
# terminated meshing, display surface
viz1 = nes_viz.DisplayDataSetSurface(faces, color=nes_msh.skin_color, display_edges=True)
nes_viz.RunRenderWindowInteractorOnDisplayedObject(viz1)
# save surface mesh as VTK format
wr = vtk.vtkPolyDataWriter()
wr.SetInputData(faces)
wr.SetFileName("GeneratedMeshTrifaces.vtk")
wr.Write()
# save volume mesh as VTK format
wr = vtk.vtkDataSetWriter()
wr.SetInputData(mesh)
wr.SetFileName("GeneratedVolumeMesh.vtk")
wr.Write()
# save mesh as matlab file
nes_msh.SaveMeshInMatlabFormat("generated_mesh.mat", mesh, faces)
# save mesh to be directory used in Eidors
nes_msh.SaveMeshInEidorsFormat("MIRAN_generated_Eidors_mesh.mat", mesh, faces)
Finally a cut-through the volume mesh is represented below using ParaView. In this representation two different colors (red and beige) have been used to represent the sub-regions corresponding to the segmented lungs. The mesh surface is represented as a wireframe in light cyan color, to provide a representation of the body volume.
ACKNOWLEDGEMENTS
MIRAN is a software library based on C++ high-performance classes developed by NE Scientific. These classes make use of the following third party software:
- VTK (Visualization Toolkit) for visualization and other purposes.
- GDCM for DICOM image loading and manipulation.
- VTK_DICOM for interfacing to GDCM from within VTK classes.
- AVCD for decimating surface meshes.
- Sparse Field Active Contours by Shawn Lankton – a modified and enhanced version of this library is used for segmentation purposes.
- TetGen is used for generating volume meshes from the segmented surfaces.