Three-dimensional numerical schemes for the segmentation of the psoas muscle in X-ray computed tomography images

Abstract

The analysis of the psoas muscle in morphological and functional imaging has proved to be an accurate approach to assess sarcopenia, i.e. a systemic loss of skeletal muscle mass and function that may be correlated to multifactorial etiological aspects. The inclusion of sarcopenia assessment into a radiological workflow would need the implementation of computational pipelines for image processing that guarantee segmentation reliability and a significant degree of automation. The present study utilizes three-dimensional numerical schemes for psoas segmentation in low-dose X-ray computed tomography images. Specifically, here we focused on the level set methodology and compared the performances of two standard approaches, a classical evolution model and a three-dimension geodesic model, with the performances of an original first-order modification of this latter one. The results of this analysis show that these gradient-based schemes guarantee reliability with respect to manual segmentation and that the first-order scheme requires a computational burden that is significantly smaller than the one needed by the second-order approach.

Introduction

Sarcopenia [1] is a pathological condition of the skeletal muscle characterized by loss of mass, and depletion of strength and physical performance. In the last decade, several studies have shown that sarcopenia is associated to specific chronic diseases. Just as examples, and not exhaustively, low skeletal muscle mass is associated with post-transplantation death [2, 3]; it is a risk factor for poor performances and toxicity of chemotherapy intervention in cancer patients [4,5,6,7,8]; it has a high probability of co-occurrence with several neurodegenerative conditions such as Parkinson’s disease and multiple sclerosis [9, 10].

Despite its relevance as a frailty biomarker, sarcopenia is difficult to quantitatively determine, most studies in this field currently relying on estimates of variation of the mass of the psoas muscle determined by post-processing specific slices of low resolution X-ray computed tomography (CT) scans [11]. Specifically, in a typical pipeline utilizing clinical data, abdominal CT images of the patient are downloaded from the hospital Picture Archiving and Communication System (PACS); a cross-sectional image in correspondence with the third lumbar vertebra (L3) is selected; the psoas is identified in this image; and the corresponding psoas area and density are computed.

Different approaches to automated two-dimensional psoas segmentation have been proposed. For example, anatomy-based manually extracted prior shapes have been utilized to initialize energy-based iterative schemes [12,13,14]. More recently, artificial intelligence (AI) relied on convolutional neural networks (CNNs) and Generative Adversarial Networks (GANs) to realize CT images segmentation specifically tailored to psoas identification at L3 [15, 16]. Of course, deep learning algorithms generally designed for muscles and internal organs segmentation [17, 18] can be applied to the specific task of psoas identification.

This image-based approach to sarcopenia quantification relies on the controversial assumption that a single CT slice is representative of the whole psoas condition [19,20,21]. Very recently [20, 22], two studies conceived in our (extended) research group accounted for this ambiguity and showed that a complete three-dimensional segmentation of the psoas muscle in CT images of oncological and neurological patients may be exploited to obtain metabolic information about the diseased tissue. Both studies utilized hybrid data made of low-dose CT and positron emission tomography (PET) data of patients injected with 18F-Fluorodeoxyglucose (FDG-PET/CT), and, in both cases, the computational paradigm was based on the following three steps [23,24,25]:

• 1. A segmentation algorithm was applied to co-registered CT data in order to automatically identify the psoas muscle.

• 2. A binary mask was constructed, where each voxel within the image region corresponding to the muscle was given value 1 and all other voxels were given value 0.

• 3. The binary mask was voxel-wise multiplied with the FDG-PET volumes in order to extract the metabolic information associated to the psoas tissue.

Using this scheme, paper [20] introduced a novel image-driven biomarker, the Attenuation Metabolic Index (AMI), that can be applied as a prognostic tool in metastatic castration-resistant prostate cancer patients, while paper [22] proved that there is a common underlying mechanism for skeletal muscle and spinal cord hypermetabolism in amyotrophic lateral sclerosis (ALS).

Of course, the reliability of these metabolic outcomes strongly depends on the accuracy with which the segmentation step 1 above is realized. Specifically, in the first study [20] the segmentation approach was based on an extended version of the Hough transform [26], which can be used only if a realistic parametric model for the muscle is at disposal; the second study [22] utilized a rather heuristic approach based on histogram equalization and an $\alpha$-shape algorithm to identify the region corresponding to the inner muscle. However, despite the increasing potential relevance of this kind of analysis for diagnostic/prognostic applications, no systematic study has been performed so far on the effectiveness of three dimensional numerical schemes for image segmentation of the psoas muscle in low-resolution CT data.

The objective of the present paper is to address the computational problem of segmenting the psoas muscle in clinical low-dose CT data by means of three dimensional numerical schemes originating from the level set methodology [27,28,29,30,31]. The use of three-dimensional methods is particularly appropriate for assessing sarcopenia. Indeed, volumetric segmentation naturally exploits the information coming from adjacent regions, improving the accuracy with which the tissue density is computed [32, 33]. Specifically, we focused on gradient-based methods that search for an optimal solution via a search orientation based on the minimization of an objective function [34]. These methods can incorporate in their objective function weighted curve length and variation of intensity in the segmented regions [35], and can be combined with deep learning techniques [36]. In the present study we considered a classical evolution model, a standard second-order geodesic model, and a first-order modification of the geodesic model that is formulated for the first time in the present paper. In all three cases the model discretization has been performed by means of a standard finite-difference scheme. As a result, we obtained three segmentation algorithms, all with the nice property that they require the users to select just one point for each psoas in the data volume to initialize the segmentation process. The application to CT images contained in hybrid PET/CT volumes showed a notable segmentation accuracy, although with rather significant differences in the required computational burden. However, due to the low spatial resolution of these clinical CT data, in some cases the segmented shapes presented some artifacts at the end of the prescribed iterations. Therefore, we have formulated and implemented an algorithm for removing these false positives, which can be applied at the end of the evolution process when necessary.

We finally point out that three-dimensional segmentation schemes allow the leveraging of information from adjacent slices [37]. However, in the case of medical imaging applications, the use of three-dimensional deep-learning-based algorithms may be hampered by the computational burden and, sometimes, by the lack of sufficiently large data sets, necessary in the training phase of the networks [32]. Therefore, the introduction of three-dimensional deterministic numerical schemes may be particularly supportive for this kind of applications.

The plan of the paper is as follows. “Evolution models and their discretization” section describes the level set models employed in the study, together with the numerical scheme utilized for their discretization. “Pre-processing and post-processing” section discusses the more heuristic approaches exploited to pre-processing the CT images and post-processing the outcomes of the numerical evolution process. “Numerical results” section compares the results provided by the different numerical procedures when applied to clinical CT data. Our conclusions are offered in “Comments and conclusions” section.

Evolution models and their discretization

The general framework considered in the present study assumes that segmentation is realized by letting a three-dimensional wave-front

$$\begin{aligned} \gamma _t:=\{x \in \mathbb {R}^3 : v(t,x)=0\} \end{aligned}$$

(1)

evolve from an initial point inside the shape to segment until it reaches its boundary. Following the level set method [27], the field $v=v(t,x)$ is assumed to be the solution of the Cauchy problem

$$\begin{aligned} \left\{ \begin{array}{ll} v_t+g(x)|D v|=0,\qquad & (x,t)\in \mathbb R^3\times \mathbb (0,T)\\ v(0,x)=v_0(x)& x\in \mathbb R^3, \end{array} \right. \end{aligned}$$

(2)

where D is the three-dimensional gradient operator, g(x) is the speed in the normal direction and $v_0$ is a proper representation of the front $\gamma _0$ (i.e. denoting by $\Omega _0$ the region enclosed by the front, $v_0$ has to be equal to 0 on $\gamma _0$, positive outside $\Omega _0$ and negative inside). It is well known that, under appropriate assumptions on the velocity field g(x) and the initial condition $v_0$, the viscosity solution of (2) exists and is unique (see, e.g., [38]). In the following we refer to the eikonal Eq. (2) as the classical model. This approach can be generalized to a second order Hamilton-Jacobi equation of the form

$$\begin{aligned} v_t+H(x,Dv,D^2v)=0, \end{aligned}$$

(3)

which includes more sophisticated definitions of the velocity field g, possibly dependent on the curvature of the surface, as for example in the Mean Curvature Motion. Here we consider the version proposed by Caselles et al. in [39] denominated geodesic model

$$\begin{aligned} H(x,Dv,D^2v) = - g(x) |D v| \left( \varepsilon \nabla \cdot \left( \frac{Dv}{|Dv|}\right) -\mu \right) - \eta Dg \cdot Dv, \end{aligned}$$

(4)

where $\varepsilon , \mu , \eta$ are positive parameters to appropriately tune. To give a better insight on the choices made in the numerical section below, we recall that these parameters are used to properly weigh the different components of the model. For instance, $\eta$ is related to the transport driven term, which attracts the front towards the boundaries of the image, and $\varepsilon$ is related to the curvature part that regularizes the solution and tends to shorten the front. Finally, the constant velocity $\mu$ is added to balance the shrinking effect of the curvature term and allows outward evolution, as in our case, where the front expands from an initial sphere towards the surface of the muscle. We decided to give a normalized value to the constant $\mu =1$ and adapted the other parameters accordingly. More precisely, the higher the value of $\eta$, the slower the expansion toward the boundary, whereas high values for the curvature term enforce the use of very restrictive Courant-Friedrichs-Lewy (CFL) conditions to keep the finite difference scheme stable.

We also propose a first-order modification of Eq. (4), i.e.

$$\begin{aligned} H(x,Dv) = \mu g(x) |D v|- \eta Dg \cdot Dv , \end{aligned}$$

(5)

which is characterized by a lower computational burden. In order to complete the description of the segmentation model we have to define the velocity g(x) such that it forces the front towards the relevant boundary and then it stops its evolution. A typical definition in gradient-based methods is the following

$$\begin{aligned} g(x)=\frac{1}{\left( 1+|\nabla (G*I(x))|^p\right) },\qquad p\ge 1, \end{aligned}$$

(6)

where I represents the intensity values of the image, and G a Gaussian filter of deviation $\sigma$, used to regularize noisy images. This function has the desired properties for edge detection; in fact, it takes values in [0, 1], and it is decreasing with respect to the gradient of the image I. If a spike in the gradient is present, then g is close to 0; otherwise it is equal to 1 in smooth regions. In practice, the filtering process is done once before the evolution, and the values of the smoothed image $G \star I$ are used instead of I in the computation of the gradient. High values for the parameter p may enlarge the region where the function is close to 0, reducing the overall accuracy of the scheme. This is why here we used (6) with $p=2$ as in [39].

We now provide some details about the numerical schemes utilized for the discretization of the first order Eqs. (2) and (5), which can also be adapted to Eq. (4) with minor changes. For sake of simplicity the description of these schemes is provided in the two-dimensional setting, although the codes we used for data analysis have been implemented for three-dimensional segmentation.

We realized the numerical discretization of all three models by means of the standard difference form

$$\begin{aligned} u^{n+1}_{i,j} = S_\Delta (u^{n})_{i,j} := u^{n}_{i,j} -\Delta t ~ h\left( x_j,y_i,D_x^-u^n_{i,j},D_x^+u^n_{i,j},D_y^-u^n_{i,j},D_y^+u^n_{i,j}\right) , \end{aligned}$$

(7)

where the numerical Hamiltonian is the local Lax-Friedrichs Hamiltonian

$$\begin{aligned} h(x,y,p^-,p^+,q^-,q^+) := & H\left( x,y,\frac{p^++p^-}{2},\frac{q^++q^-}{2}\right) \nonumber \\ & -\frac{\alpha _x(p^-,p^+)}{2}(p^+-p^-)-\frac{\alpha _y(q^-,q^+)}{2}(q^+-q^-), \end{aligned}$$

(8)

with

$$\begin{aligned} \alpha _x(p^-,p^+) := \underset{\underset{p\in I(p^-,p^+)}{x,y,q,}}{\textrm{max}}\left| H_p(x,y,p,q)\right| ,\qquad \alpha _y(q^-,q^+) := \underset{\underset{q\in I(q^-,q^+)}{x,y,p,}}{\textrm{max}}\left| H_q(x,y,p,q)\right| , \end{aligned}$$

(9)

with $I(a,b) := [\min (a,b),\max (a,b)]$. We point out that $h(x,y,p^-,p^+,q^-,q^+)$ is a Lipschitz continuous function, with

$$\begin{aligned} D_x^{\pm } u^n_{i,j}:=\pm \frac{u^n_{i,j\pm 1}-u^n_{i,j}}{\Delta x}, \end{aligned}$$

(10)

and

$$\begin{aligned} D_y^{\pm }u^n_{i,j}:=\pm \frac{u^n_{i\pm 1,j}-u^n_{i,j}}{\Delta y}. \end{aligned}$$

(11)

Further, in this context consistency is imposed via condition

$$\begin{aligned} h(x,a,a,b,b)=H(a,a,b), \end{aligned}$$

(12)

whereas monotonicity requires the numerical Hamiltonian h to be non-decreasing with respect to its second and fourth argument and non-increasing with respect to the third and the fifth ones. The resulting scheme (7)-(8) is monotone under the Courant-Friedrichs-Lewy (CFL) type condition $\frac{\Delta t}{\Delta x}\cdot \alpha _x + \frac{\Delta t}{\Delta y}\cdot \alpha _y \le 1$. Moreover, as suggested in [40], this scheme can also be used to solve Eq. (4) by first splitting the Hamiltonian into two components

$$\begin{aligned} H(x,y,Dv,D^2 v)=H_0(x,y,Dv)+H_1(k(x,y),Dv), \end{aligned}$$

where $H_0$ includes the eikonal and advection term, whereas $H_1$ accounts for the dependence on the curvature k(x, y), which is computed as

$$\begin{aligned} k(x,y)=\nabla \cdot \left( \frac{Dv(x,y)}{|Dv(x,y)|}\right) = \frac{v_{xx} v_y^2 - 2v_xv_yv_{xy}+v_{yy}v_x^2}{{(v_x^2+v_y^2)}^{3/2}}. \end{aligned}$$

Then, it is enough to discretize the first-order part using the Lax-Friedrichs scheme and $H_1$ with simple centered finite difference approximations. In this case, the stability of the scheme requires a parabolic-type CFL condition of the form $\Delta t \sim \Delta x^2$.

Pre-processing and post-processing

Before solving the evolution equations, we pre-processed the CT data according to the following process:

• 1. Hounsfield normalization: given that normal bone tissue in CT has a Hounsfield unit (HU) of around 160, all HU values above 120 have been set to 0, in order to reduce the relative weight of the bones [41].

• 2. For image enhancement purposes, a histogram equalization algorithm has been applied, which performed a re-mapping of the image based on the probability distribution of the input gray levels [42].

• 3. A Gaussian filter with standard deviation $\sigma$ has been applied to improve the signal-to-noise ratio of each image [43].

This process was executed slice-wise in a two dimensional manner. Finally, we created a mask for all HU values inside the interval [5, 120] to save the position of the non-muscular tissue, which is needed to improve the stability of the scheme for practical application.

The initialization of the scheme requires the user to select two points, one for each muscle. These points become the centers of the two spheres of radius $r=5\Delta x$ used as initial condition for the evolution, which is executed separately for each psoas. For all simulations the spatial discretization (the length of the voxel side) is set to $\Delta x = 0.1$, whereas the maximum number of iterations $N_{max}$ is heuristically chosen case by case, depending on the number of slices in the CT volume and the overall quality of the images.

As previously remarked, clinical low-dose CT images as the ones acquired by PET/CT scanners, are usually noisy and present a lower level of resolution with respect to high-resolution CT acquisitions. Further, the complexity of the human body in the psoas region often implies the presence of sections where different organs or muscles, which have close HU values, have numerous contact points. Therefore, the definition of a clear boundary of the psoas becomes an intricate issue, which impacts the reliability of the numerical scheme in automatically segmenting the image. Specifically, in some cases it may happen that the wave-front leaves the correct boundary and starts detecting spurious objects outside the muscle of interest, thus making the usual stopping criterion ineffective, with unreliable final results.

In order to solve this issue, we devised a post-processing routine able to clean the result of (most of) the external objects and eventually give an acceptable final result. The procedure is composed of two steps and has a very low computational cost, which makes it efficient even if more than one iteration of the second step is needed. The routine works separately for each slice, but makes use of the information given by the closest slice in order to exploit the intrinsic continuity of the surrounding problem. We have summarized the procedure in Algorithm 1, while more details on the implementation are given in the following items.

• First step. The main aim of this step is to reduce the number of connected components (CCs) in each slice, detecting the one belonging to the psoas and discarding the remaining ones. To do so, the algorithm starts from the slice containing the center chosen by the user and selects the correct CC by trivially excluding the ones that do not include the chosen point. Then, the remaining slices are scanned toward the top and the bottom, separately, always starting from the center. The main loop consists in first computing the centroid of the CC detected in the previous slice, and then, in the current slice, in selecting just the CC containing such point. The final result is a segmented muscle in which only the spurious objects directly connected to the boundary may be still present.

• Second step. In this second step the starting slice and the previous one, in the direction of the scanning, are assumed not to need the cleaning procedure. With this assumption, a vector tol is initialized with length equal to P, the set of segmented voxels in the starting slice, where the $i-$th entry is the sum of minimum distance of the voxel $i \in P$ from the set of voxels $P_{prev}$ in the previous slice with a fixed constant (0.75). Then the algorithm loops over the slice, setting $P_{prev} = P$ and first computing the difference set $P_{diff}=P-P_{prev}$, where P and $P_{prev}$ are the sets of voxels in the current and previous slices, respectively; then, for each point $j \in P_{diff}$, the algorithm computes the set distance $dist(j,P_{prev})$ with respect to $P_{prev}$. This distance is obtained in correspondence of a voxel $i_{min} \in P_{prev}$. If $dist(j, P_{prev})>tol(i_{min})$, voxel j is excluded from P and hence removed from the segmentation. After all voxels in $P_{diff}$ have been scanned and removed if needed, a new tolerance vector tol is computed using the current cleaned set of voxels P and $P_{prev}$.

The second step, in its initialization phase, assumes that the starting slice and the previous one (according to the scanning direction) do not need any cleaning. If this is not the case, a second iteration of this step is needed, using a different starting point.

We finally point out that, on the one hand, the pre-processing procedures described at the beginning of this section, which are rather standard, can be realized by means of algorithms different than the ones employed in this study. For example, image enhancement can be obtained by means of methods based on integral transforms of AI-based techniques (see [44] for a review of these methods); and noise reduction can be obtained by means of adaptive or anisotropic Gaussian filtering [45, 46]. On the other hand, the post-processing procedure has been tailored on the three-dimensional schemes introduced in this study, although both steps address issues that are typical of image segmentation performed by means of level set methods. At a more general level, we validated the outcomes of the overall segmentation process with respect to modifications of the input parameters in the pre-processing steps, and found out a notable robustness of the outcomes. However, for a few subjects such outcomes (before post-processing) are characterized by a lower quality (see Fig. 3 below for some examples), which is due to the topographical complexity of the district close to the psoas and to the poor resolution of the low-dose CT volumes. For such cases, the second step of the post-processing phase typically requires more than one iteration. In any case, we think that a more systematic investigation of the impact of different post-processing procedures on the segmentation quality is worthwhile consideration for future studies.

Algorithm 1 Cleaning procedure

Numerical results

The pipeline made of pre-processing, numerical segmentation, and post-processing has been applied to the CT images (see Fig. 1) acquired for nine subjects belonging to a retrospective data set of patients recruited at the IRCCS Ospedale Policlinico San Martino, Genova, Italy. All subjects were submitted to FDG-PET/CT and the study was performed according to the Declaration of Helsinki, Good Clinical Practice, and local ethics regulation. All enrolled patients signed a written informed consent at the time of FDG-PET/CT, econmpassing the use of anonymized data for retrospective research purposes.

For numerical segmentation we used the three numerical schemes introduced in “Evolution models and their discretization” section, i.e., the finite difference scheme (7)-(11) applied to (2) (‘classical’ from now on), and to the two geodesic schemes (4) and (5) (‘GMFD 1ord’ and ‘GMFD 2ord’ from now on, respectively). For both the two finite difference goedesic model schemes we set $\mu =1$ and $\nu =0.25$, while $\epsilon =0.05$ for the second order model. Further, ’GMFD 1ord’ exploits the linear CFL condition $\Delta t = \lambda \Delta x$, with

$$\begin{aligned} \lambda = (3 \Vert g \Vert + \Vert D_x g \Vert + \Vert D_y g \Vert + \Vert D_z g\Vert )^{-1}, \end{aligned}$$

(13)

and ‘GMFD 2ord’ utilizes the parabolic CFL condition $\Delta t = 0.25 \Delta x^2$.

Fig. 1

Coronal view of the psoas muscles (in yellow) and the corresponding axial views at different levels of the abdomen

Table 1 contains the number of iterations and CPU time employed to realize the segmentation by the three numerical schemes when applied to the volumes of all nine patients. Figure 2 contains the results of the segmentation provided by the three methods in the case of six subjects and Fig. 3 describes the performances of the post-processing algorithm in the case of the remaining three subjects. Specifically, the three cases in this latter figure are the ones for which the outcomes of the numerical schemes present a particularly significant amount of artifacts. However, the post-processing step can be applied to all outcomes, including the ones described in Fig. 2. The segmentations in Figs. 2 and 3 show that the run of the numerical schemes without any post-processing may imply either the persistence of some artifacts that typically show up in the region between the two psoas muscles, or a slight shrinking of the psoas volume. The presence of false positive or shrinking may impact the assessment of sarcopenia realized by means of the PET/CT-based process described in the Introduction. In fact, an inaccurate evaluation of the psoas dimension implies the construction of inaccurate binary maps that are used to extract metabolic information from the PET data, and may result in an over- or under-estimated determination of the FDG uptake of the psoas tissue.

Table 1 Number of iterations, CPU time, values of p and $\sigma$, and number of slices for the three segmentation algorithms for the nine CT volumes considered in this study

We compared the segmentation results we obtained with the three methods described in the previous sections with the ones provided by manually drawn profiles within the OsiriX software package [47] and by the 3D Slicer extension TotalSegmentator (TS) [48], a fully automated total body segmentation deep learning routine based on nnU-Net [49], which has been applied with reliable results for the segmentation of abdominal tissues even in the three-dimensional setting [50,51,52]. In our analysis, TS has been used in both fast (TS-FM) and full-resolution (TS-FRM) modes. In order to evaluate the agreement between the manually drawn psoas profiles and the segmentations obtained with the three numerical schemes and the TotalSegmentator tool we used six metrics: the Dice similarity coefficient [53], the Jaccard index [54], the Hausdorff distance [55], the Average Symmetric Surface Distance (ASSD) [56], the Surface Dice (Surface Dice) and the Robust Hausdorff distance (Robust Hausdorff). Surface Dice and Robust Hausdorff have been computed via the packages provided by https://github.com/google-deepmind/surface-distance/tree/master.

The Dice similarity coefficient and the Jaccard index vary between 0 and 1 where 0 represents an empty intersection between the two segmentations while 1 a perfect overlap. The Hausdorff distance and ASSD are indeed distances, so smaller values correspond to closer objects. Specifically, the Hausdorff distance is the greatest distance from a point in one segmentation mask to the closest point in the other mask, while ASSD measures the average of all distances from points on the boundary of the first mask to the boundary of the second one, and viceversa. Similarly to the classical volumetric Dice, Surface Dice measures the overlap of two surfaces, where a surface element is counted as overlapping when the closest distance to the other surface is less than or equal to the specified tolerance, fixed at 3 mm. The range of the Surface Dice is [0,1]. Robust Hausdorff is computed like the classical Hausdorff distance but it uses a fixed percentile of the distances instead of the maximum distance. We choose the 75th percentile. In Table 2 we show the average values of these metrics computed across the nine subjects, together with the corresponding standard deviations.

Fig. 2

Segmentations provided by the classical evolution model (2) (‘classical’, first column), the first-order geodesic model (5) (‘GMFD 1ord’, second column), and the second order geodesic model (4) (‘GMFD 2ord’, third column). First row: subject CTTO02 (45 slices) with $p=2$ and $\sigma =0$. Second row: subject PZGE01 (33 slices) with $p=2$ and $\sigma =2$. Third row: subject PZGE02 (35 slices) with $p=2$ and $\sigma =2$. Fourth row: subject PZGE03 (31 slices) with $p=2$ and $\sigma =2$. Fifth row: subject PZTO03 (52 slices) with $p=2$ and $\sigma =0$. Sixth row: subject PZTO04 (51 slices) with $p=2$ and $\sigma =1$

Fig. 3

Post-processing pipeline for the patients CTTO01, CTTO03, and PZTO01 when applied to the outcomes of ‘GMFD 1ord’. First column: the starting data. Second colum: the outcome of the first step (starting slice: 20 for CTTO01, 44 for CTTO03, 25 for PZTO01). Third column: the outcome of the second step

Table 2 Reliability assessment for the three three-dimensional numerical schemes when compared to manual segmentation. The metrics utilized for this analysis are the Dice coefficient (optimal value: 1), the Surface Dice coefficient (optimal value: 1), the Jaccard coefficient (optimal value: 1), the Hausdorff distance, the Robust Hausdorff distance, and ASSD distance. The last two columns report the metrics values in the case of two standard software tools utilized in medical imaging, i.e., the TotalSegmentator 3D extension of Slicer in both fast and full-resolution mode (TS-FM and TS-FRM, respectively)

We also tested the above described techniques on data downloaded from the Cancer Imaging Archive at https://www.cancerimagingarchive.net/. Specifically, we used six subjects collected in the Cancer Moonshot Biobank - Prostate Cancer Collection (CMB-PCA) [57, 58]. Figure 4 shows the segmentations of a dicom image of the psoas for one of the six subject in the dataset, while the average values of the metrics computed across the six subjects are shown in Table 3.

We finally point out that, on the one hand, the computational burden required by the ‘classical’ and ‘GMFD 1ord’ schemes is one order of magnitude smaller than the one required by ‘GMFD 2ord’. On the other hand, deep learning segmentation in fast mode (TS-FM) is typically performed within a computational time that is almost one order of magnitude smaller than the ones employed by the two fastest level set schemes (although with a rather significantly lower segmentation reliability).

Fig. 4

Different segmentations for one of the six subjects in the CMB-PCA dataset. First row, from left to right: ground truth, TS-FM and TS-FRM segmentations. Second row, from left to right: ‘classical’, ‘GMFD 1ord’ and ‘GMFD 2ord’ segmentations

Table 3 Reliability assessment for the three three-dimensional numerical schemes when compared to manual segmentation on the six subjects in the CMB-PCA dataset. The metrics utilized for this analysis are the Dice coefficient (optimal value: 1), the Surface Dice coefficient (optimal value: 1), the Jaccard coefficient (optimal value: 1), the Hausdorff distance, the Robust Hausdorff distance, and ASSD distance. The last two columns report the metrics values in the case of two standard software tools utilized in medical imaging, i.e., the TotalSegmentator 3D extension of Slicer in both fast and full-resolution mode (TS-FM and TS-FRM, respectively)

Comments and conclusions

This paper has proposed a 3D approach to the segmentation of the psoas muscle based on level set models. These numerical schemes provide a higher degree of automation with respect to the approach utilized in [22], which is notably heuristic, and a higher degree of generality with respect to the approach utilized in [20], whose performances rely on the availability of prior geometrical models for the psoas. For all numerical tests the geodesic models provide stable results, whereas the classical scheme may produce over-resolved shape, as it is particularly evident in the cases of patient CTTO02 (Fig. 2, first column, first row) and patient PZGE02 (Fig. 2, first column, third row). This effect is not evident for the two geodesic models, probably because the transport term in (4)-(5) pushes the front toward the external boundary of the muscle. Further, this same transport term usually helps the front stick to the right boundary, reducing the overflow outside the muscle, which happens when some complex structure is present in some portion of the data, as shown in Fig. 2, fifth row. As far as the reliability of the segmentation is concerned, the average Hausdorff and ASSD distances computed between the masks produced by our schemes and the manually identified ground-truth are, respectively, more than two times and more than 1.5 times lower than the distances obtained using the nnUNet-based TS masks (see Tables 2 and 3). This fact highlights that, even if the Jaccard and the Dice metrics indicate comparable overlap (of both volume and surface), our segmentation masks present less outlier. These results are confirmed by the use of the Surface Dice modification, which is $4-5$% higher with respect to TS, and of the Robust Haussdorff distance, which is 22% lower, on average. In the case of the images downloaded from the CMB-PCA dataset, TS and the three level set numerical schemes present a similar degree of reliability if Dice, Surface Dice and Jaccard indicators are applied, while the deterministic numerical methods have a significantly better behaviour than deep learning according to the Haussdorff and ASSD metrics (see Table 3).

The first order version of the geodesic model has performances with comparable reliability to the ones characterizing the second order model. Indeed, the metrics and distances used in Table 2 have similar mean values and uncertainties. These same uncertainties are clearly smaller than the ones associated to the classical model. Further, the computational cost required by ’GMFD 1ord’ to achieve its metrics and distances values is one order of magnitude smaller than the one typical of the second order model (see Table 1). The higher computational burden of ’GMFD 2ord’ is mainly due to the parabolic CFL condition, which slows down the evolution and forces the algorithm to utilize a greater number of iterations to reach convergence.

It is clear that the simple post-processing procedure devised to clean the results provided by the numerical schemes is able to remove all spurious objects in most cases, still keeping the volume of the segmented muscle mostly untouched with respect to the results of segmentation. From a heuristic viewpoint, we could test that, if the outcome given by the scheme is already accurate on its own, one iteration of the procedure is enough to remove most artifacts. On the other hand, if the spurious objects are many and, in particular, some of them intersect the slice containing the center of the initial condition, a second run with a (possibly) different starting point is necessary. In any case, if the latter point is chosen correctly, the final result can be satisfactory and only in few cases some small objects still persist.

From a computational viewpoint, we observe that higher order schemes, as defined for example in [59], could be considered and adapted to the modified model (4) to increase the accuracy of the segmentation at the expense of longer computational times, but this goes beyond the scope of the present work and will be considered in future research.

The clinical impact of the present study can be of potential interest. In fact, in general, a systematic application of three-dimensional approaches for the segmentation of psoas can contribute to assess the limitations of the currently widespread clinical protocol, which is based on the analysis of a single cross-sectional image in correspondence to L3. More specifically, the use of deterministic numerical scheme like the ones based on level set methods may have some advantage with respect to AI-based models, since machine/deep learning may suffer data scarcity in the training phase.

In conclusion, on the one hand, this approach to segmentation of medical imaging data can be of course extended to the case of other muscular districts, although the psoas dimension and shape make this muscle particularly appropriate for automated segmentation in the case of sarcopenia assessment. On the other hand, the application of this approach in the case of other imaging modalities like MRI and PET can be problematic, due to the lower inter-organ contrast typical of MRI data and low spatial resolution typical of PET data. Finally, future work could include exploring adaptive parametric selection in the numerical schemes, integration with deep learning models, and applications to prospective clinical studies.