Introduction

Measuring the thickness of an anatomical structure is a common task for radiologists. More specifically, in the field of cardiology, it is often needed to measure the thickness of the myocardial muscle, as it convey information about its health status, which has been shown to be related to its electrophysiological properties (Ghannam Michael et al. (2018),Takigawa et al. (2019)).

The visible myocardium segmented on medical images can have inclusions of non-muscle tissue such as fat infiltration, calcifications, or fibrosis. It can happen that we are interested in measuring the thickness of the actual muscle fibers, ignoring these inclusions (fig. 1). In this article we describe a method to perform this task, and show that it can reciproqually be used to measure the thickness of the inclusion.

Methodology

Our method is an extension of the method described by Yezzi and Prince (Yezzi and Prince (2003)).

Correspondence trajectories

Briefly, Yezzi and Prince define thickness at each point \(x\) of the tissue region \(R\) as the total arclength of a unique curve, passing through \(x\). This curve originates on the the inner boundary of the tissue region \(\partial_0 R\), and terminates on its outer boundary \(\partial_1 R\) (see fig. 1)

Such curves can be obtained by solving the Laplace equation over R:

\[ \Delta u = 0 \](1)

with the Dirichlet boundary conditions:

\[ u(\partial_0 R) = 0 \textnormal{~and~} u(\partial_1 R) = 1\](2)

\(u\) defines a scalar field in \(R\) and is used to define curves within \(R\) with the tangent field \(\overrightarrow{T}\):

\[ \overrightarrow{T} = \frac{\nabla u}{|| \nabla u ||} \](3)

For didactic purposes, in fig. 1 the curve is represented as a straight dotted line (“thickness”), but \(\overrightarrow{T}\) actually defines curvy trajectories for complex shapes of \(R\).

Arclength computation

Yezzi and Prince (Yezzi and Prince (2003)) define the length functions \(L_0\), where \(L_0(x)\) gives the arclength of the correspondence trajectory between \(\partial_0 R\) and \(x\) (reciprocally, \(L_1\)(x) between \(\partial_1 R\) and \(x\)).

\[ \nabla L_0 \cdot \overrightarrow{T} = 1, \textnormal{~with~} L_0(\partial_0 R) = 0 \](4)

\[ -\nabla L_1 \cdot \overrightarrow{T} = 1, \textnormal{~with~} L_1(\partial_1 R) = 0 \](5)

Thickness at every point \(x\) is then defined by summing these length functions:

\[ W(x) = L_0(x) + L_1(x) \](6)

Weighted arclength

To account for the inclusion of tissue which thickness we want to exclude, we propose to replace the right term (constant) of eq. 4 and eq. 5 by a function \(f\) varying over the tissue region.

\[ \nabla L_0 \cdot \overrightarrow{T} = f(x) \](7)

\[ -\nabla L_1 \cdot \overrightarrow{T} = f(x) \](8)

We call \(f\) the arclength weight function and it has the following properties:

• If x is fully occupied by the tissue region we want to measure, \(f(x) = 1\)
• If x is fully occupied by an inclusion we want to exclude from the thickness measure, \(f(x) = 0\)

\(f(x)\) can also take any value in the \([0, 1]\) interval. This can be useful to account for partial volume effects by using a transfer function between the intensity (eg, Hounsfield Units) of the medical image and \(f\).

Ray-tracing as a possible alternative

Another approach to achieve comparable measurements is to use ray-tracing. In such framework, one needs to define straight lines emanating from one surface and measure the length of the segment between the two surfaces \(\partial_0 R\) and \(\partial_1 R\). Typically, these lines would be normal to the surface of the considered anatomical structure at the point where they emanate. The measured distance (ie, the tissue thickness) could be similarly weighted with the weight function \(f\).

However, our methodology presents several advantages over this approach:

• It operates directly on voxels and it does not require to define a surface, something usually achieved by converting voxel data to triangular meshes. While there are several algorithms and tools to achieve this conversion, arbitrary choices must be made to control, among other things, the smoothness of the surface; and these choices, influencing the directions of the rays, can have drastic effects on the thickness values.
• The results with our method are thickness values for all voxels of the segmented structure, and not on a single surface. This can be leveraged for further analysis, or parameterizing fast electrophysiological models based on rectangular-grid data (Rapaka et al. (2012), Cedilnik et al. (2017)).
• There is no need to define out of which surface these trajectories have to be calculated. Since normals have different directions depending on whether they are defined from the inner or outer surface, this would result in visually different thickness maps depending on the considered surface. With our method, and as a consequence of both the definition of the correspondence trajectories and the voxel nature of the output, one can see topologically similar thicknesses on the epicardial, endocardial, or even mid-wall surface.

Implementation

Algorithm

It has been shown (Yezzi and Prince (2003)) that computing eq. 4 and eq. 5 over \(R\), in 3D amounts to solving (equations (8) and (9) of Yezzi and Prince (2003))

\[ L_0[i, j, k]=\frac{f(x)+\left|T_x\right| L_0[i \mp 1, j, k]+\left|T_y\right| L_0[i, j, \mp 1, k]+\left|T_z\right| L_0[i, j, k \mp 1]}{\left|T_x\right|+\left|T_y\right|+\left|T_z\right|} \]

\[ L_1[i, j, k]=\frac{f(x)+\left|T_x\right| L_1[i \pm 1, j, k]+\left|T_y\right| L_1[i, j, \pm 1, k]+\left|T_z\right| L_1[i, j, k \pm 1]}{\left|T_x\right|+\left|T_y\right|+\left|T_z\right|} \]

NB: in Yezzi and Prince (2003), \(f(x)\) is a constant equal to 1.

Numerical solving

Results can be obtained by the same methods described by Yezzi et Prince (Yezzi and Prince (2003)). Iterative approaches until convergence are possible, but a fast-marching-like algorithm is more efficient (computational resource-wise), especially in a single-threaded context, since it requires a single pass over \(R\), in theory. However we want to note here, than unlike Yezzi and Prince (2003) describe, our experiments showed that this fast-marching-like approach does not reach convergence in a single pass. We propose an iterative/fast-marchine-like combined approach where the traversal order of the first pass is memorized for subsequent passes, reaching convergence more rapidly than a naive iterative approach.¹

Results on toy geometries

One dimension

We conducted experiments in one dimension to verify that the results were those expected. We set correspondence trajectories to straight horizontal lines and set a 4-element wide “hole” along these trajectories where the weight function values ranged from 0 to 1 (fig. 2).

As expected, the resulting weighted thickness \(W_p(x)\) equals to the number of elements \(n\) in R along the horizontal line when the weight function \(f(x)\) is set to 1 along the line, corresponding to the non-weighted thickness computation. When \(f(x_{\textnormal{holes}}) = 0\), \(W_p(x) = n - n_{\textnormal{holes}}\). For \(0 < f(x_{\textnormal{holes}}) < 1\), \(n < W_p(x) < - n_{\textnormal{holes}}\).

Two dimensions

For illustrative purposes, before showing applications on clinical data in 3D, we also present a 2 dimensional toy geometry (fig. 3). It resembles a radiological “short-axis view” of a human heart, in which the myocardium presents inclusion of pixels where \(f(x) < 1\). As can be seen on fig. 3 (rightmost picture), the parts of the (toy) myocardium on the left where such inclusions are present present a smaller weighted thickness values compared to the surrounding, “non-weighted” regions.

Applications

In this section we present possible applications of our method on clinical data.

Workflow

For all clinical applications, the preliminary step is a segmentation of the structure that we want to measure. For the myocardium, this is typically done by segmenting the blood pool inside the cardiac cavities to create an endocardial binary mask, and the outer layer of the myocardial wall to create an epicardial mask (Komatsu et al. (2013)).

One then has to define the weight function \(f\). The way it is defined varies depending on the application, as shown in the following sections.

Muscle-only thickness

Fat and calcification inclusions in the myocardium are common. Evaluating the thickness of the muscle fibers excluding these inclusions is the most straightforward application of our method.

The Houndsfield units (HU) of muscle and fat are close, especially in the context of iodine-injected CT. After segmentation, we performed a histogram analysis of HU values within the myocardium to determine the peak value \(P\), and standard deviation \(\delta\). We define a transfer function to set \(f(x)\) depending on the HU such that:

• for \(I(x) > P - \delta\), \(f(x) = 1\)
• for \(I(x) < P - 2\delta\), \(f(x) = 0\)
• for \(2\delta < I(x) < \delta\), \(f(x)\) is linearly interpolated

An example result is shown on fig. 4.

Scar thickness

Our method can also be used to quantify the local thickness of an infarct scar. To achieve this, after segmentation of the endocardial and epicardial masks, one has to segment the scar mask (how to obtain such segmentation is outside the scope of this article).

The weight function \(f\) is then defined such that \(f(x) = 0\) in the myocardium and \(f(x) = 1\) in the scar. This can be used to visually assess the thickness of the fibrosis, as shown on fig. 5

Scar transmurality

Instead of the absolute thickness of the scar, in millimeters, our method can also be used to visualise the local proportion of fibrosis in the myocardial wall, a ratio we call scar transmurality. A value of 1 would mean that an area of the myocardial wall is pure fibrosis and reciproqually, a value of 0 must be interpreted as the absence of fibrosis in this area, ie. no ischemic scar.

To obtain such measure, a possible approach is to:

1. Compute the non-weighted thickness of the wall, ie, setting \(f(x) = 1\) all over the myocardial wall.
2. Perform the voxel-wise ratio of scar thickness (see previous section) over the total thickness.

The transmural extent, i.e., the scar transmurality is defined for every voxel of the myocardial wall. See fig. 6 for an example result.

Conclusion

In this article, we defined a measure we named weighted thickness, suited to measure the thickness of a tissue on medical images, when it presents inclusions of a different tissue that we want to exclude from the measure. It does so by giving voxels different weights depending on a weight function that must be defined accordingly to the considered application. Like other partial differential equation-based thickness measures, it outputs a scalar map of the same dimensions as the input segmentation, in which each pixel (or voxel) has a thickness value, coherent from the inner to the outer boundary of the measured structure.

The main limitation of our article is that we do not evaluate the clinical relevance of the weighted thickness measure. In particular, the transfer function described in \(\ref{sec:muscle}\) has been arbitratily defined and would require tuning and validation for the implied application, ie defining the health status of the myocardial on CT images. Our future work will focus on such validation.