Current framework is organized in four interconnected modules summarized in Fig. 1: Pre-processing, Depression quantification, Inner chest contour segmentation and Thoracic indexes computation. The software code has been developed in MATLAB® 2020a (https://it.mathworks.com/), running in Windows 10.
As a preliminary step for subsequent analyses, a range of slices of interest must be selected, including the slice of maximal sternal depression, on which measurements for PE indices are usually performed in clinic. Indeed, not all axial images acquired are useful for our analysis, but only the ones in which the chest depression and the lungs are clearly visible and thus could be quantified. Of course, this excludes marginal slices at the beginning and at the end of the scan, where the amount of depression is negligible.
This step has been implemented through a Graphical User Interface for selecting range of slices, slice for PE indexes computation, as well as patient’s sex (Additional file 1: Fig. 1).
Pre-processing
In order to improve low contrast inherent to MR images, firstly we perform a contrast adjustment by remapping the values of the input intensity to fill the entire intensity range. Then, we focus exclusively on chest district by excluding arms placed at the borders of images, due to the small dimension of chest in pediatric patients. This is obtained by defining a proper mask, based on subject’s thorax morphology (Additional file 1: Fig. 2).
Depression quantification
This module has the goal to quantify the depression, based on a volumetric study. Indeed, rather than evaluating the depression on a single slice, as traditional radiological indices commonly adopted in clinical practice do, we propose to analyze multiple slices in order to measure the depression volume. The idea is to identify the two maximum and the minimum points of the outer chest contour for each slice considered and thus define an elliptic curve between the two maximum points to correct the depression and simulate the normal chest, in absence of PE malformation. The difference between the chest image before and after image correction gives the amount of the depression.
Analysis of outer chest contour
First of all, the algorithm turns the grey-scale image into a binary image, by applying a manual threshold (T = 0.1) to separate the foreground from the background pixels. Indeed, this value proves to apply for all examined subjects. Then, we get exterior boundaries of chest in terms of Cartesian x and y coordinates through morphological gradient operators in order to identify the two maximum and the minimum points of outer chest contour. This is just a rough detection of maximum points, since binarization may alter upper profile of images. We thus resort to morphological operators of closing to correct upper image boundaries. By only considering the upper half of outer chest boundary, the algorithm exploits the spatial discontinuities of boundary pixels along y direction in order to identify the two maximum points (Fig. 2a,b). As regards the minimum point, we use boundary pixel locations before morphological operations and find it as lowest peak in the range of y coordinates between the two already identified maximum points (Fig. 2c–e).
Depression volume
By analyzing MR images of chest in normal patients, we have noticed that the best curve, representing the chest morphology between the two maximum points of outer chest contour, could be a roto-translated ellipse, whose points are calculated as follow:
$$\begin{aligned} & x = \frac{{x_{1} + x_{2} }}{2} + a*cos\left( t \right)*cos\left( \alpha \right) - b*sin\left( t \right)*sin\left( \alpha \right) \\ & y = \frac{{y_{1} + y_{2} }}{2} + a*cos\left( t \right)*sin\left( \alpha \right) - b*sin\left( t \right)*cos\left( \alpha \right) \\ \end{aligned}$$
where:
-
(x, y): coordinates of ellipse points
-
(x1, y1): coordinates of right vertex of major axis
-
(x2, y2): coordinates of left vertex of major axis
-
a: semi-major axis
-
b: semi-minor axis, found as \(b=a*\sqrt{1-{ e}^{2}}\)
-
e: ellipse eccentricity
-
t: variation angle of ellipse points, defined between 0 and π (half ellipse)
-
α: rotation angle, defined as the angle between the horizontal line and major axis
Eccentricity (e) and positions of right and left vertices of major axis ((x1, y1) and (x2, y2)) represent the parameters we have modified in order to simulate the profile of outer chest contour in normal patients. After several tests, e has been set to 0.99.
Regarding the position of major axis vertices, we have separated patients based on their sex. Indeed, anatomical differences between male and female forced us to deal with the depression issue in a distinct way. By analyzing normal chest images of male subjects, we were able to find a unique method to define the position of major axis vertices. Specifically, the algorithm identifies y coordinates (y1 and y2) by lowering the position of two maximum points of outer chest contour by a constant value, while the x coordinates (× 1 and × 2) are found by searching the most extreme points at the same y coordinates (Fig. 3a).
This correction method does not work in female patients, because the presence of breast rises the position of major axis vertices by causing a wrong depression correction. As among female subjects there is a high variability in chest shape due both to age and differentiated anatomical growth, it is impossible to define a single correction method for the depression. For these reasons, we decided to exclude female patients from depression quantification analysis, while including them in the automatic computation of standard indexes.
After the ellipse has been obtained, the algorithm finds the indices corresponding to x and y ellipse coordinates in the image matrix and adds pixels to the binary image of chest in these specific locations. Consequently, the depression is filled by applying a morphological operation of closing using a disk as structural element (Fig. 3b).
The depression area is thus calculated as the difference between the image before application of morphological operators and the one after correction with the elliptical curve. Finally, the depression volume is computed by summing the volumes obtained for each slice, resulting from the product of depression areas by ‘slice thickness’ image DICOM attribute. After repeating this operation for each slice, the algorithm is able to represent on grey-scale images how the normal morphology of chest should be (Fig. 3c).
New pathological marker computation
The absolute value of depression volume cannot be used as a pathological marker since it is strongly dependent on chest dimension and on the number of slices considered for its computation, which is different from patient to patient. Thus, we decided to normalize it on the thorax volume after the correction, as it simulates the ‘normal’ condition of the chest. Specifically, the algorithm quantifies the correct chest volume in the same way as for depression volume by considering the binary image after depression correction. The new pathological marker, that we named Volumetric Correction Index, is defined as follow:
$$Volumetric\,Correction\,Index= \frac{depression\,volume}{correct\,chest\,volume} * 100$$
Therefore, the new index proposed represents the percentage of depression that must be corrected in PE patients.
Inner chest contour segmentation
This module aims at detecting the inner contour of the chest, fundamental for PE indexes calculation. If this task is difficult in CT images, where the attenuation coefficients of the heart and the chest-wall are quite close, it is even more challenging in MR images, where different chest regions often have a high similarity in terms of grey levels.
Therefore, in order to simplify the segmentation process, we designed this module by subdividing it in consecutive steps, as described in the following sections.
Firstly, the algorithm isolates the inner chest portion by exploiting lung segmentation and similarity between the inner and outer wall contour. Then, it excludes the vertebral body by thresholding method. Finally, it corrects the errors in the detection of inner chest contours through a comparison among consecutive slices. For this analysis step, we opted for working on a limited number of slices, by excluding the ones preceding the slice selected for index computation. Indeed, the remaining range of slices ensures an easily implementable segmentation of inner chest region thanks to an optimal lung-background contrast.
Lung segmentation
In view of performing segmentation of the lungs, we used histogram analysis for identification of the correct threshold. The grey-level histogram of a MR image is characterized by a high variability, both across subjects and across slices within the same patients, in peaks’ shapes corresponding to the lungs and to cardiac structures and thorax tissue, respectively. For this reason, Otsu thresholding technique [34], the standard approach for histogram partitioning, does not perform well due to its inability to correctly separate bimodal histograms when the two classes are very different in size. Therefore, we developed a method to automatically partition a grey-level histogram, by adapting an algorithm presented by [35]. The idea proposed by this study, that we applied to our problem, is to locate the concavity between the two principal peaks in the curve representing the image histogram by maximizing divergence between the histogram and a Gaussian fit.
After computing the histogram of the grey-scale image in continuous form, the algorithm defines an auxiliary curve P(x) on the same grey-level range of the histogram H(x). We assumed P(x) as a normal distribution, with mean given by µ, the average gray-level of H(x), and the corresponding variance given by σ2. We also considered P(x) and H(x) to have an identical area α under their curves. Given that xmin ≤ x ≤ xmax, P(x) is defined as:\(P(x) = \frac{\alpha }{z} G(x)\)where:
-
\(G(x)= \frac{1}{\sqrt{2\pi {\sigma }^{2}}}exp \left(\frac{{-(x - \mu )}^{2}}{2{\sigma }^{2}}\right)\)
-
\(z (x) =\sum_{x={x}_{min}}^{{x}_{max}}G (x)\)
-
\(\alpha =\sum_{x={x}_{min}}^{{x}_{max}}H (x)\)
-
\(\mu =\frac{1}{\alpha }\sum_{x={x}_{min}}^{x={x}_{max}} x H(x)\)
-
\({\sigma }^{2} =\frac{1}{\alpha }\sum_{x={x}_{min}}^{x={x}_{max}} {(x - \mu )}^{2}H(x)\)
As a normal distribution, P(x) presents its largest value at x = µ and has a convex part that goes from µ − σ to µ + σ. The oddity is that the highest peaks of H(x) are close to µ, the average grey level, such that the concavities surrounded by the highest peaks in H(x) are often in contrast with the convex part of P(x). Hence, the line (l) that divides the main peaks in H(x) can be easily found by maximizing the difference between P(x) and H(x), or formally: \(l~ = arg\mathop {\max }\limits_{x} ~\left( {{\text{P}}\left( x \right){\text{~}}{-}{\text{~H}}\left( {\text{x}} \right)} \right){\text{~}}\), with \(\mu - \sigma \le x\le \mu +\sigma\).
The l value corresponds to the threshold separating the two main peaks in image histogram. We implemented this histogram partitioning method in MATLAB and applied it twice in our analysis. First, it is used to separate the chest area from the background. Thus, it analyzes the lower part of histogram, by considering as input the range of x values in the low grey-level region. Once found the threshold (lbg) that removes the background from the image, a new grey-level histogram H’(x) is generated, considering only the pixels related to the chest. Hence, the method is reapplied to the new data to estimate the correct threshold for the lung segmentation (llung). Specifically, in order to enhance the threshold search, the algorithm focuses the analysis on the lower part of H’(x), since it corresponds to grey values belonging to lungs (Additional file 1: Fig. 3).
After finding appropriate thresholds for each slice with this strategy, these are used to segment lungs from chest region. Specifically, a mask is created where pixels with intensity above the llung are set as white and the remaining ones are set as black. Other segmented elements besides the lungs are removed and morphological operation of closing are applied to smooth edges and fill holes inside the lungs.
Selection of appropriate slice for indices computation
Finally, before proceeding to inner contour detection, we created an automatic technique to exclusively select the slices where lungs are clearly visible. Indeed, we wanted to exclude from further analysis those slices where inner chest region segmentation could be complex, due to the absence of lungs. Specifically, the algorithm computes the lung area and relates it to the entire chest area. Then, it selects only the slices in which the ratio is greater than 20%. Therefore, if the slice selected for indices calculation shows high similarity in grey values between different chest areas, the algorithm automatically picks the first following slice, where inner chest contour detection can be performed properly.
Inner chest contour detection
In order to isolate the inner thoracic region, we adapted an algorithm proposed by [36], for the inner curvature detection of CT images. Specifically, they proposed a recursive algorithm that exploits outer wall contour as a starter point for inner contour segmentation, due to similarity in morphology between them.
We identified as algorithm inputs, obtained from previous module, the matrices composed of pixel locations of each lung and the matrix containing pixel locations of outer curvature.
The algorithm goes through steps along the outer curvature in clockwise direction until the start point is found again. Every 12 steps the actual point and the point 12 steps before are connected and a perpendicular line in the mid-point of their connection is generated. Then the algorithm finds the intersection point between the perpendicular line and the first point crossed by it on the two lungs. In the area of the binary image where the perpendicular lines do not cross any lungs, a correction of the invalid points generated is necessary. Thus, it calculates for each line the distances between the midpoint and the intersection point and computes their mean value (µd) and the standard deviation (σd). We have designed a length filter by defining as invalid the intersection points whose distances from the mid-points are longer than a specific threshold that we identified as 2* µd − σd. All points, corresponding to longer distances than this value, are deleted and replaced by new ones located at the same distances as the previous valid point (Fig. 4a). Once all intersection points have been found, the inner curvature is calculated by an interpolating process. Initially, the algorithm prepares intersection points by separating them in two subsets: the ones related to the upper half of the inner contour and those belonging to the lower part. Additionally, it performs an initial correction, by deleting points whose y locations are in discontinuity with y positions of neighboring points, in order to avoid the eventual errors made by previous operations. Then it applies a shape-preserving piecewise cubic interpolation method (‘pchip’) with a high sampling rate. Finally, we obtain a group of closely spaced points both for superior and inferior half of inner contour. After re-combining them in a unique set of points, we can define the boundary of a mask that isolates the inner thoracic region (Fig. 4b). However, this mask also includes the vertebral body and is not accurate in all the slices, mostly due to bad lung segmentation. For these reasons, it is necessary to improve the inner chest segmentation with further processing.
Inner chest contour correction
The first step of inner chest contour correction consists in excluding the vertebral body. For doing so, the algorithm applies a thresholding method by using the inner mask found in previous section as a tool to improve segmentation. After masking out the external chest area, a slice-wise threshold is defined to exclude out heart and other cardiac elements through histogram partitioning method presented in previous section. Indeed, correction is not possible without masking out cardiac district, due to high similarity in grey values between heart and thoracic tissue. Thus, the algorithm returns to original chest image, before the application of inner mask, and assigns to background the just segmented pixels belonging to cardiac structures (Fig. 5a). Then, it is able to separate the inner chest region from the outer chest one, by applying as threshold the same value found for lung segmentation (llung). To have the inner region as foreground, the complementary image is computed, and some morphological operations are applied to smooth the edges and fill the holes. Finally, we obtain a binary image representing the inner chest region from which to extract boundary pixel locations for each slice.
However, a further correction may be necessary both around the vertebral body, having grey values close both to lungs and to cardiac structure intensities, and around inferior lung area, being often difficult separating pixels belonging to lungs and to thoracic tissue ones (Fig. 5b).
We thus designed a method to correct the inner chest contour by comparing consecutive slices thanks to high similarity in boundary pixel positions belonging to the lower half of inner contour of our interest. At this step, user intervention is required such that correction process starts from a slice where inner contour detection does not present errors. Once first slice is selected, the algorithm starts the pair-wise comparison among adjacent slices in both directions, by taking as reference the points belonging to the contour of slice selected. Thus, the algorithm computes the distance among them and all the points belonging to the contour of the adjacent slice, which is assumed as incorrect. Then it creates a vector that, for each point belonging to the contour to correct, only keeps the minimum distance among all those just computed. It also calculates the maximum value (dmax) and the standard deviation (σd) of all minimum distances. After several tests, we established that the algorithm must continue only if σd is greater than 1.8. Additionally, we identified dmax -2* σd as threshold value that separates correct points from incorrect ones. Thus, for each incorrect point, the algorithm finds the range that must be deleted, by identifying its extremities in the nearest points to the correct curve. Then, it replaces them with the points belonging to the correct curve by using a shape-preserving piecewise cubic interpolation method (‘pchip’) (Fig. 5c). Once a new curve is obtained, the algorithm proceeds to the next slice, by taking as reference the just corrected contour. Such an algorithm allows to satisfactorily correct errors in segmentations (Fig. 5d).
Thoracic indexes computation
This module aims at computing PE indices used by physicians to classify the severity of patients’ malformation. As mentioned above, among multiple thoracic markers, we focused on the severity (Haller index and Correction index) (Fig. 6a, b) and deformity (Asymmetry index and Flatness index) ones (Fig. 6c, d). The algorithm only works on the first slice of images processed in the previous module. Indeed, it corresponds to the slice selected by the user or to the first following one where inner chest contour can be detected.
Once inner distances and thoracic indices are computed, the framework saves their results along with the new pathological marker obtained in Depression quantification module in an Excel® file, located in the same folder as patient’s images. Each quantified distance in following computations has been multiplied by ‘pixel spacing’ attribute to have measures in mm.
Haller index
Haller index (iHaller) is calculated by dividing the transverse diameter, i.e., the widest horizontal distance of the inside of the ribcage, to the minimum anteroposterior diameter (min APd), i.e., the shorter distance between the vertebral body and the sternum [12].
$$iHaller=\frac{transverse\,diameter}{min\,APd}$$
As regards transverse diameter, the algorithm identifies its first extremity as the point on inner chest contour with minimum x coordinate and the second one as the point at its same y coordinate. Conversely, the first extremity for measuring min APd corresponds to sternum position. We approximate it as the point with maximum y coordinate (y values decrease toward the bottom of image) by only considering the range of inner chest contour points between x coordinate positions of two maximum points of outer chest contour upper half. Second extremity corresponds to the vertebral body position. It is taken as the point with minimum y coordinate by only considering the range of inner chest contour points between position of x coordinates of two maximum points of outer chest contour lower half.
Correction index
Correction index (iCorrection) is calculated by dividing the amount of defect, measured as the difference between the maximum anteroposterior distance, i.e. the maximum distance between the anterior spine and the anterior portion of the chest (max APd) and the minimum anteroposterior diameter (min APd), to the maximum anteroposterior distance (max APd), multiplied by 100 [37].
$$iCorrection=\frac{(max\,APd - min\,APd)}{min\,APd} * 100$$
For max APd computation, firstly, the algorithm draws a horizontal line at the same y coordinate of vertebral body position, that is assumed as the anterior spine position. Then, it identifies two points on inner chest contour at the same x coordinate of two maximum points of outer chest contour. The latter are assumed as the positions of right and left anterior portion of the chest. Thus, for each point it computes the distances between them and the horizontal line and gets the maximum diameter between the two distances.
Asymmetry index
Asymmetry index (iAsymmetry) is calculated by dividing the longest anteroposterior distance of the right chest wall (right hemithorax APd) to the longest anteroposterior distance of the left chest wall (left hemithorax APd), multiplied by 100 [38].
$$iAsymmetry=\frac{right\,hemithorax\,APd}{left\,hemithorax\,APd} * 100$$
Right hemithorax APd’s extremities are identified as the points on inner chest contour at the same x coordinate of first maximum point, as it is located in right hemithorax. Right hemithorax APd’s bounds are identified as the points on inner chest contour at the same x coordinate of second maximum point, as it is situated in left hemithorax.
Flatness index
Flatness index (iFlatness) is computed by dividing the transverse diameter of the thorax to the longer of the two maximum anteroposterior diameters of the right (right hemithorax APd) and left hemithorax (left hemithorax APd) [13]. As all the distances have been already found, the algorithm can proceed with Flatness index computation, as follow:
$$iFlatness=\frac{transverse\,diameter}{max\left(righthemithoraxAPd;\,lefthemithoraxAPd\right)}$$
User’s correction
As mentioned above, the algorithm does not always perform indices computation on the same slice selected by user, due to its inability to segment images where different chest areas have similar grey values. In these cases, it selects the first following slice, where inner chest contour detection can be performed. We noticed that by going through consecutive slices some inner distances maintain their value constant, while others, specifically min APd and max APd are more likely to vary. For this reason, in the algorithm we add the possibility of user’s intervention when the slice is different from the one selected. Specifically, the user is asked to insert two points on the image, useful for min APd and max APd calculation: sternum position and vertebral body position. Thus, the algorithm recomputes the indices by considering the modifications on these two inner distances. Finally, two sets of results are obtained: the ones calculated on the slice picked by the algorithm and those obtained on the same slice selected after correction of sternum and vertebral body points.