Advertisement

A patient-specific approach for quantitative and automatic analysis of computed tomography images in lung disease: Application to COVID-19 patients

  • L. Berta
    Affiliations
    Department of Medical Physics, ASST Grande Ospedale Metropolitano Niguarda, Piazza Ospedale Maggiore 3, Milan 20162, Italy
    Search for articles by this author
  • C. De Mattia
    Affiliations
    Department of Medical Physics, ASST Grande Ospedale Metropolitano Niguarda, Piazza Ospedale Maggiore 3, Milan 20162, Italy
    Search for articles by this author
  • F. Rizzetto
    Affiliations
    Department of Radiology, ASST Grande Ospedale Metropolitano Niguarda, Piazza Ospedale Maggiore 3, Milan 20162, Italy
    Search for articles by this author
  • S. Carrazza
    Affiliations
    Department of Physics, Università degli Studi di Milano and INFN Sezione di Milano, via Giovanni Celoria 16, Milan 20133, Italy
    Search for articles by this author
  • P.E. Colombo
    Affiliations
    Department of Medical Physics, ASST Grande Ospedale Metropolitano Niguarda, Piazza Ospedale Maggiore 3, Milan 20162, Italy

    Department of Physics, Università degli Studi di Milano and INFN Sezione di Milano, via Giovanni Celoria 16, Milan 20133, Italy
    Search for articles by this author
  • R. Fumagalli
    Affiliations
    Department of Medicine and Surgery, University of Milan-Bicocca, Monza, Italy

    Department of Anaesthesia and Intensive Care Medicine, ASST Grande Ospedale Metropolitano Niguarda, Piazza Ospedale Maggiore 3, Milan 20162, Italy
    Search for articles by this author
  • T. Langer
    Affiliations
    Department of Medicine and Surgery, University of Milan-Bicocca, Monza, Italy

    Department of Anaesthesia and Intensive Care Medicine, ASST Grande Ospedale Metropolitano Niguarda, Piazza Ospedale Maggiore 3, Milan 20162, Italy
    Search for articles by this author
  • D. Lizio
    Affiliations
    Department of Medical Physics, ASST Grande Ospedale Metropolitano Niguarda, Piazza Ospedale Maggiore 3, Milan 20162, Italy
    Search for articles by this author
  • A. Vanzulli
    Affiliations
    Department of Oncology and Hemato-Oncology, Università degli Studi di Milano, via Festa del Perdono 7, Milan 20122, Italy

    Department of Radiology, ASST Grande Ospedale Metropolitano Niguarda, Piazza Ospedale Maggiore 3, Milan 20162, Italy
    Search for articles by this author
  • A. Torresin
    Correspondence
    Corresponding author.
    Affiliations
    Department of Medical Physics, ASST Grande Ospedale Metropolitano Niguarda, Piazza Ospedale Maggiore 3, Milan 20162, Italy

    Department of Physics, Università degli Studi di Milano and INFN Sezione di Milano, via Giovanni Celoria 16, Milan 20133, Italy
    Search for articles by this author
Published:January 28, 2021DOI:https://doi.org/10.1016/j.ejmp.2021.01.004

      Highlights

      • A new approach has been proposed to overcome limitations of quantitative metrics calculated from lung CT images.
      • New metrics, derived from physics assumptions and with physiological significance, are introduced.
      • New metrics show lower dependencies from CT-number and reduced inter and intra patient variability.
      • Quantitative metrics of lung COVID-19 diseases can be described by subdividing the organ into 24 sub-regions.

      Abstract

      Purpose

      Quantitative metrics in lung computed tomography (CT) images have been widely used, often without a clear connection with physiology. This work proposes a patient-independent model for the estimation of well-aerated volume of lungs in CT images (WAVE).

      Methods

      A Gaussian fit, with mean (Mu.f) and width (Sigma.f) values, was applied to the lower CT histogram data points of the lung to provide the estimation of the well-aerated lung volume (WAVE.f). Independence from CT reconstruction parameters and respiratory cycle was analysed using healthy lung CT images and 4DCT acquisitions. The Gaussian metrics and first order radiomic features calculated for a third cohort of COVID-19 patients were compared with those relative to healthy lungs. Each lung was further segmented in 24 subregions and a new biomarker derived from Gaussian fit parameter Mu.f was proposed to represent the local density changes.

      Results

      WAVE.f resulted independent from the respiratory motion in 80% of the cases. Differences of 1%, 2% and up to 14% resulted comparing a moderate iterative strength and FBP algorithm, 1 and 3 mm of slice thickness and different reconstruction kernel. Healthy subjects were significantly different from COVID-19 patients for all the metrics calculated. Graphical representation of the local biomarker provides spatial and quantitative information in a single 2D picture.

      Conclusions

      Unlike other metrics based on fixed histogram thresholds, this model is able to consider the inter- and intra-subject variability. In addition, it defines a local biomarker to quantify the severity of the disease, independently of the observer.

      Keywords

      1. Introduction

      COVID-19 is a complex infectious disease characterized by common and non-specific symptoms, such as fever, cough, shortness of breath and fatigue, and a broad spectrum of clinical manifestation, ranging from asymptomatic infection to respiratory failure requiring oxygen support or invasive ventilation [
      • Zhou F.
      • Yu T.
      • Du R.
      • Fan G.
      • Liu Y.
      • Liu Z.
      • et al.
      Clinical course and risk factors for mortality of adult inpatients with COVID-19 in Wuhan, China: a retrospective cohort study.
      ].
      Computed tomography (CT) is the current standard of reference to assess lung alteration, even at early stage of the disease, when the patient has few or no symptoms, and to monitor the course of the disease at different time points [
      • Chung M.
      • Bernheim A.
      • Mei X.
      • Zhang N.
      • Huang M.
      • Zeng X.
      • et al.
      CT imaging features of 2019 novel coronavirus (2019-NCoV).
      ,
      • Rubin G.D.
      • Ryerson C.J.
      • Haramati L.B.
      • Sverzellati N.
      • Kanne J.P.
      • Raoof S.
      • et al.
      The Role of Chest Imaging in Patient Management During the COVID-19 Pandemic.
      ]. Despite the increase in chest CT examinations due to COVID-19 pandemic, the use of low dose protocols guarantees a very low risk of cancer induction [
      • Ghetti C.
      • Ortenzia O.
      • Maddalo M.
      • Altabella L.
      • Sverzellati N.
      Dosimetric and radiation cancer risk evaluation of high resolution thorax CT during COVID-19 outbreak.
      ].
      The majority of COVID-19 studies use a qualitative approach, describing the lesions by visual and pictorial assessment. A lexicon for the description of chest CT imaging findings in coronavirus disease, i.e. the COVID-19 Reporting And Data System (COVID-RADS) [
      • Salehi S.
      • Abedi A.
      • Balakrishnan S.
      • Gholamrezanezhad A.
      Coronavirus disease 2019 (COVID-19) imaging reporting and data system (COVID-RADS) and common lexicon: a proposal based on the imaging data of 37 studies.
      ], was proposed in addition to the terminology endorsed by the Fleischner Society Nomenclature Committee [
      • Hansell D.M.
      • Bankier A.A.
      • MacMahon H.
      • McLoud T.C.
      • Müller N.L.
      • Remy J.
      Fleischner Society: Glossary of terms for thoracic imaging.
      ,

      S. Simpson F.U. Kay S. Abbara S. Bhalla J.H. Chung M. Chung et al. Radiological Society of North America Expert Consensus Statement on Reporting Chest CT Findings Related to COVID-19. Endorsed by the Society of Thoracic Radiology, the American College of Radiology, and RSNA. Radiol Cardiothorac Imaging 2020;2:e200152. 10.1148/ryct.2020200152.

      ]. The purpose is to standardize terminology and communication and assessing a possible COVID-19 presence.
      In addition, visual quantitative analysis, i.e. the scoring of lung abnormalities assessed by visual interpretation of CT images, and the densitometric evaluations based on the histogram of the Hounsfield Unit (HU) distribution, were demonstrated useful to predict clinical severity [
      • Li K.
      • Fang Y.
      • Li W.
      • Pan C.
      • Qin P.
      • Zhong Y.
      • et al.
      CT image visual quantitative evaluation and clinical classification of coronavirus disease (COVID-19).
      ,
      • Yang R.
      • Li X.
      • Liu H.
      • Zhen Y.
      • Zhang X.
      • Xiong Q.
      • et al.
      Chest CT Severity Score: An Imaging Tool for Assessing Severe COVID-19.
      ,
      • Sumikawa H.
      • Johkoh T.
      • Yamamoto S.
      • Yanagawa M.
      • Inoue A.
      • Honda O.
      • et al.
      Computed tomography values calculation and volume histogram analysis for various computed tomographic patterns of diffuse lung diseases.
      ].
      Lung segmentation, the first step required for a quantitative analysis of medical images [

      M. Mascalchi G. Camiciottoli S. Diciotti Lung densitometry: Why, how and when J Thorac Dis 9 2017 3319 45 https://doi.org/10.21037/jtd.2017.08.17.

      ], can be performed nowadays through several fully automatic tools that help to reduce the intra and the inter-reader variability [
      • Wang J.
      • Li F.
      • Li Q.
      Automated segmentation of lungs with severe interstitial lung disease in CT.
      ,
      • Park B.
      • Park H.
      • Lee S.M.
      • Seo J.B.
      • Kim N.
      Lung Segmentation on HRCT and Volumetric CT for Diffuse Interstitial Lung Disease Using Deep Convolutional Neural Networks.
      ].
      Many indexes can be derived from the lung density histogram, starting from simple measurements, as the mean density value [
      • Nakagawa H.
      • Nagatani Y.
      • Takahashi M.
      • Ogawa E.
      • Van Tho N.
      • Ryujin Y.
      • et al.
      Quantitative CT analysis of honeycombing area in idiopathic pulmonary fibrosis: Correlations with pulmonary function tests.
      ], to the measurement of the relative area of emphysema in patients with chronic obstructive pulmonary diseases [

      Lynch D. Progress in Imaging COPD, 2004 - 2014. Chronic Obstr Pulm Dis J COPD Found 2014. https://doi.org/10.15326/jcopdf.1.1.2014.0125.

      ] or the extraction of descriptive parameters of the histogram such as kurtosis and skewness [
      • Ash S.Y.
      • Harmouche R.
      • Vallejo D.L.L.
      • Villalba J.A.
      • Ostridge K.
      • Gunville R.
      • et al.
      Densitometric and local histogram based analysis of computed tomography images in patients with idiopathic pulmonary fibrosis.
      ].
      Analysis of CT images has already been used in the past to better understand the pathophysiology of the acute respiratory distress syndrome (ARDS) [
      • Gattinoni L.
      • Caironi P.
      • Pelosi P.
      • Goodman L.R.
      What Has Computed Tomography Taught Us about the Acute Respiratory Distress Syndrome?.
      ], using the density histogram to define lung compartments with different aeration.
      In literature, different values of Hounsfield Units (HU) have been proposed to define these compartments [
      • Gattinoni L.
      • Pesenti A.
      • Bombino M.
      • Baglioni S.
      • Rivolta M.
      • Rossi F.
      • et al.
      Relationships Between Lung Computed Tomographic Density, Gas Exchange, and PEEP in Acute Respiratory Failure.
      ]. To quantify the well-ventilated regions in COVID-19 patients, a potential surrogate to estimate the residual respiratory function and the alveolar recruitment during ventilation [
      • Nishiyama A.
      • Kawata N.
      • Yokota H.
      • Sugiura T.
      • Matsumura Y.
      • Higashide T.
      • et al.
      A predictive factor for patients with acute respiratory distress syndrome: CT lung volumetry of the well-aerated region as an automated method.
      ,
      • Gattinoni L.
      • Caironi P.
      • Cressoni M.
      • Chiumello D.
      • Ranieri V.M.
      • Quintel M.
      • et al.
      Lung recruitment in patients with the acute respiratory distress syndrome.
      ], the interval between −950 HU and −700 HU was proposed [
      • Chen A.
      • Karwoski R.A.
      • Gierada D.S.
      • Bartholmai B.J.
      • Koo C.W.
      Quantitative CT analysis of diffuse lung disease.
      ,

      D. Colombi F.C. Bodini M. Petrini G. Maffi N. Morelli G. Milanese et al. Well-aerated Lung on Admitting Chest CT to Predict Adverse Outcome in COVID-19 Pneumonia Radiology 2020:201433. 10.1148/radiol.2020201433.

      ]. This threshold was already used in the past for studies involving other lung diseases [
      • Matsuoka S.
      • Yamashiro T.
      • Matsushita S.
      • Kotoku A.
      • Fujikawa A.
      • Yagihashi K.
      • et al.
      Quantitative CT evaluation in patients with combined pulmonary fibrosis and emphysema: Correlation with pulmonary function.
      ,
      • Ohkubo H.
      • Kanemitsu Y.
      • Uemura T.
      • Takakuwa O.
      • Takemura M.
      • Maeno K.
      • et al.
      Normal lung quantification in usual interstitial pneumonia pattern: The impact of threshold-based volumetric CT analysis for the staging of idiopathic pulmonary fibrosis.
      ], even if without a fully consensus, especially for ARDS. For example, a recent work used the density of −500 HU, already suggested by Gattinoni et al. [
      • Gattinoni L.
      • Pesenti A.
      • Bombino M.
      • Baglioni S.
      • Rivolta M.
      • Rossi F.
      • et al.
      Relationships Between Lung Computed Tomographic Density, Gas Exchange, and PEEP in Acute Respiratory Failure.
      ], to discriminate between normally (−900, −501 HU) and compromised lung (-500, +100 HU), then further divided in poorly (−101, −500 HU) and non-aerated (−100, +100 HU) [

      Lanza E, Muglia R, Bolengo I, Santonocito OG. Quantitative Chest CT analysis in COVID-19 to predict the need for oxygenation support and intubation. PrePrint 2020:1–18. https://doi.org/10.21203/rs.3.rs-30481/v1.

      ].
      Despite being extremely informative, chest quantitative analysis has several drawbacks limiting its routine application.
      First, the natural inter-patient and intra-patient variability, mainly due to the respiratory cycle. Several studies performed inspiratory and/or expiratory breath holds in combination with mechanical ventilation to standardize lung inflation during image acquisition [
      • Gattinoni L.
      • Caironi P.
      • Cressoni M.
      • Chiumello D.
      • Ranieri V.M.
      • Quintel M.
      • et al.
      Lung recruitment in patients with the acute respiratory distress syndrome.
      ,
      • Chiumello D.
      • Langer T.
      • Vecchi V.
      • Luoni S.
      • Colombo A.
      • Brioni M.
      • et al.
      Low-dose chest computed tomography for quantitative and visual anatomical analysis in patients with acute respiratory distress syndrome.
      ]. In clinical practice, however, chest scans are obtained at the breathing-in point, inviting the patient to hold the breath. Holding the breath may however be a struggle for patients with a compromised status [
      • Cressoni M.
      • Gallazzi E.
      • Chiurazzi C.
      • Marino A.
      • Brioni M.
      • Menga F.
      • et al.
      Limits of normality of quantitative thoracic CT analysis.
      ] and it is difficult to monitor respiratory phase without dedicated equipment.
      Another limitation of current quantitative approaches is that the results are often not readily interpretable as physiological parameters but only as mathematical descriptors of the distribution of voxel values of a digital image.
      Finally, acquisition and reconstruction CT parameters inevitably affect any kind of quantitative imaging analysis, with an impact on its usefulness which has to be evaluated for each individual source of variability [
      • Traverso A.
      • Wee L.
      • Dekker A.
      • Gillies R.
      Repeatability and reproducibility of radiomic features: a systematic review.
      ].The use of a different reconstruction kernel than Visual CT (VCT) studies and the standardization of the protocols for the studies dedicated to Quantitative CT analysis (QCT) of the lungs was recommended by Neweel et al.[
      • Newell Jr., J.D.
      • Tschirren J.
      • Peterson S.
      • Beinlich M.
      • Sieren J.
      Quantitative CT of interstitial lung disease.
      ].
      The aim of this study was therefore to develop an automatic patient-customized lung analyzer to overcome the definition of aerated and pathologic regions based on thresholding of the density histogram, taking into account the patient variability and different CT acquisition protocols. Different reconstructions of the same scan series of healthy lung CT images, changing algorithm, kernel and the slice thickness, were used to test the reproducibility and the uncertainty of the proposed method. A dataset of 4DCT images, used routinely for guiding radiotherapy treatment planning [
      • Tian Y.
      • Miao J.
      • Liu Z.
      • Huang P.
      • Wang W.
      • Wang X.
      • et al.
      Availability of a simplified lung ventilation imaging algorithm based on four-dimensional computed tomography.
      ], allowed to better understand the impact of the respiratory cycle on the density histogram and the related metrics. Finally, the analysis was applied on a COVID-19 cohort to distinguish the well-aerated and the pathological regions. A comparison with a method of CT lung histogram analysis based on pre-determined thresholds was added.

      2. Material and methods

      This retrospective study was approved by the Local Ethics Committee. The need for informed consent from individual patients was waived owing to the retrospective nature of the study.

      2.1 Study population

      We analyzed three cohorts of patients.
      The first one included 20 patients (10 males, 10 females, median age 47.9 years, range 14–93 years), with a CT scan including the entire lung volume but not for pulmonary diseases (healthy lungs without alterations). The exams were retrospectively collected in ASST Niguarda Hospital in July-August 2020 from Emergency Department CT scanner.
      The second cohort was composed of 20 4DCT of locally-advanced non-small cell lung cancer patients (tumor volume: median 51 cc, range 7–392 cc), publicly available in the Cancer Imaging Archive [
      • Clark K.
      • Vendt B.
      • Smith K.
      • Freymann J.
      • Kirby J.
      • Koppel P.
      • et al.
      The Cancer Imaging Archive (TCIA): maintaining and operating a public information repository.
      ], acquired at the VCU Massey Cancer Center in the Department of Radiation Oncology, from 2008 through 2012, as a reference image for radiotherapy planning.
      The third cohort was composed of 20 patients (17 males, 3 females, median age 58.5 years, range 33–73 years) randomly selected amongst those admitted in intensive care unit within the 48 h after CT scan acquisition in March 2020 in Niguarda Hospital, with positive CT chest and positive Real-Time Polymerase Chain Reaction for SARS-CoV-2.

      2.2 CT protocol

      2.2.1 Dataset 1

      CT studies of the first cohort were all acquired on a single CT scanner, a Somatom Edge unit, (Siemens AG, Forchheim, Germany) and with the same acquisition protocol. CT scans involving the entire lung were performed using a whole-body protocol with contrast agent and automatic exposure control (AEC) and automatic selection of the tube voltage (14 studies at 120 kVp, 5 at 100 kVp, 1 at 140 kVp). We selected the basal CT phase, contrast not-enhanced, reconstructed for diagnostic intent (VCT series), using iterative algorithms (IR, Safire, S1), 3 mm as slice thickness and sharp kernel (Bl57).

      2.2.2 Dataset 1a

      Using the raw data of the first patient cohort, several sets of reconstructed images were obtained changing the slice thickness (1, 3, and 5 mm), the kernel (Bl57, Br38) and the strength of the iterative algorithm from pure Filtered Back Projection (FBP) to different level of SAFIRE blending (IR-S1, IR-S3, IR-S5).

      2.2.3 Dataset 2

      The second cohort 4DCT images were acquired on the CT 16-slice Brilliance Big Bore (Philips Medical Systems, Andover, MA), in helical mode at 120 kVp, acquiring the respiratory signal trough the Real-time Position Management of the Varian Medical Systems. The raw data was sorting in 10 breathing phases, identified as a percentage, where the 0% phase corresponds to end of inhalation. Each 3D image was reconstructed with 3 mm slice thickness range using a soft kernel. All information regarding 4DCT protocol and patients of this database can be found in reference [
      • Hugo G.D.
      • Weiss E.
      • Sleeman W.C.
      • Balik S.
      • Keall P.J.
      • Lu J.
      • et al.
      A longitudinal four-dimensional computed tomography and cone beam computed tomography dataset for image-guided radiation therapy research in lung cancer.
      ].

      2.2.4 Dataset 3

      For the third cohort, a CT protocol on the same scanner of dataset 1, with the same acquisition and reconstruction parameters, was used. Non-contrast chest CT scans were performed in supine position, during inspiratory breath-hold, within the limits of the collaboration status of the patient.
      The scan parameters for each dataset are detailed in Table 1. Table 2 summarizes the CT reconstruction series added to the included exams.
      Table 1Demographic description of the datasets used and acquisition parameters of the CT series. For each dataset, the investigated variables were reported.
      DatasetStudyNumber of patients (M/F)Age (y) [Range]Scanner ModelVendorTube Voltage (kV)Parameters Studied
      1Emergency access without pulmonary disease20 (10/10)48.5[14–93]Somatom EdgeSiemens120 (n=14), 100 (n=5), 140 (n=1)Fit parameters
      1aEmergency access without pulmonary disease20(10/10)48.5[14–93]Somatom EdgeSiemens120Kernel, slice thickness and algorithm
      24DCT20(ND)NDBrilliance Big MorePhilips120Respiratory phase
      3Emergency access with RT-PCR positive20 (17/3)58.5[33–73]Somatom EdgeSiemens120Biomarkers for patients with compromised lung
      M: Male; F: Female; RT-PCR: Real-Time Polymerase Chain Reaction
      Table 2Details of the reconstructed series calculated from the raw data of the dataset 1.
      DatasetStudyVendorSlice Thickness (mm)Reconstruction AlgorithmConvolution Kernel
      1aEmergency access without pulmonary diseaseSiemens1Safire S1Bl57
      3Safire S1Bl57
      5Safire S1Bl57
      3FBPBl57
      3Safire S3Bl57
      3Safire S5Bl57
      3Safire S1Br38

      2.3 Images analysis

      2.3.1 Lung segmentation

      Anonymized datasets were exported to a dedicated workstation where, through the extension module Chest Imaging Platform (CIP, Applied Chest Imaging Laboratory; Boston, Massachusetts, USA) of the open-source 3D Slicer software (version 4.10.2, https://www.slicer.org) [

      A. Fedorov R. Beichel J. Kalpathy-Cramer J. Finet J.-C.C. Fillion-Robin S. Pujol et al. 3D Slicer as an image computing platform for the Quantitative Imaging Network Magn Reson Imaging 2012;30:1323–41. 10.1016/j.mri.2012.05.001.

      ], a fully automatic segmentation of the lung was performed. For each patient we segmented the first reconstruction series and then we applied the same mask at different slice thickness and kernel. For the dataset 2, each 3D series representing a phase of the 4D acquisition was segmented.
      Lung segmentations of COVID-19 dataset images were carefully reviewed by an experienced radiologist and manually corrected where automatic algorithm failed.
      CIP module automatically distinguishes the right and the left lung, and provides a subdivision of lung in upper, middle and inferior regions, using anatomical landmarks. As a first step, lung analysis was performed considering the whole pulmonary volume.
      Subsequently, the imaging features of the COVID-19 dataset were analyzed and interpreted in detail through the further segmentation of lung regions in subregions. We added the ventral-dorsal and medial-distal subdivision, for a total of 24 subregions.
      In each axial slice of the upper and middle lung region, the line connecting the centroids of the right and of the left lung was used to separate ventral and dorsal regions. The distal and medial regions were identified by the perpendicular lines passing through the centroids of each lung [Fig. 1]. For the inferior region, the subdivision was calculated extending the results of the lowest slice in the middle lung region.
      Figure thumbnail gr1
      Fig. 1Graphical representation of the lung volume subdivided into: left, right, superior, middle, inferior, ventral, dorsal, medial, distal, for a total of 24 sub-volumes.

      2.3.2 Well-Aerated volume estimation

      A dedicated software was written in JavaScript language to automatically calculate and analyze the histograms of the CT images within ImageJ environment [

      Rasband W.S. ImageJ U.S. National Institutes of Health 2018 Bethesda, Maryland, USA. Http://ImagejNihGov/Ij/.

      ]. This software needs as input a 3DCT chest acquisition and applies the mask of the lung and its sub volumes. All the following analysis concerns only the voxels identified by the lung mask. A relative frequency histogram of HU values is then calculated in the range −1020 ÷ +300 HU with bin width chosen by the user.
      The proposed method estimates the well-aerated volume (WAVE) under the assumption that the distribution of the HU values of the voxels regarding exclusively healthy parenchymal tissue would be described by a Gaussian function [
      • Kemerink G.J.
      • Kruize H.H.
      • Lamers R.J.S.
      • van Engelshoven J.M.A.
      Density resolution in quantitative computed tomography of foam and lung.
      ]. In fact, the healthy parenchymal tissue can be seen as a mix of air and water arranged in cellular architecture with an average uniform density. Furthermore, due to the random nature of cellular architecture, when measured at a voxel scale, the density is no longer constant. Indeed, healthy lung densitometry has a distribution of different values and the resulting density histogram would be similar to a Gaussian distribution function. If the main source of noise in CT images is due to quantum noise, the expected distribution of HU values in cellular material is still described by a Gaussian curve.
      The lung segmentation encompasses mainly the healthy parenchyma, but other structures are included as well. For instance, vessels and/or lesions are frequently included in the region of interest. These structures, due to their higher physical density, extend the tail of the histogram towards higher values of HU. However, if the altered tissue covers a limited volume, the left side of the histogram could still be representative of healthy lung tissue only and have a Gaussian shape centered on the histogram modal value.
      Based on these assumptions, WAVE.f is defined by integrating between −1020 to + 300 HU the Gaussian function with formula:
      yx=Height.fe-x-Mu.f22Sigma.f2
      (1)


      used to fit the points around the modal value in a CT lung histogram.
      Mu.f and Sigma.f are the central and the standard deviation values of the Gaussian curve, respectively, and Height.f is a normalization factor, not relevant in this work.
      Two approaches were used to obtain the fit parameters using the Curve Fittter Tool implemented in ImageJv.1.53a setting 6000 as maximum iterations number, 2 as number of restarts and 10-10 as error tolerance. A polynomial fit of the second order was used to fit the logarithm of the HU histogram frequencies with the formula:
      ln(y)=a+bx+cx2
      (2)


      Alternatively, the fit parameters were calculated with the using the “Gaussian (no offset) function”, readily available between the fitting functions.
      The relationships between parameters of Eqs. (1) and (2) are the following:
      Heigt.f=e(a-b24c)


      Mu.f=-b2c


      Sigma.f=-12c


      The fitting range was defined as the HU points in the histogram included within a starting point, on the left of the modal value, and an end point, on the right of the modal value, identified by their relative frequency values (Fig. 2).
      Figure thumbnail gr2
      Fig. 2Examples of histogram analysis. Data are normalized to the modal value which is indicated with a vertical gray line. The black line represents the total lung histogram; the blue lines show the extremes of the range of the data points used for fitting, reported as yellow crosses. Example A), CT lung histogram of a COVID-19 patient. Example B) CT lung histogram of a patient with emphysema from dataset 2. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
      We set: 3, 5, and 10 as bin width in the histogram; 20%, 30% and 40% as percentage of the modal value for the starting points and 80%, 85% and 90% of the modal value for the end point. A total of 27 combinations of bin width and range values were explored for the Gaussian fitting of the CT lung histogram points on the dataset 1 in order to evaluate the WAVE.f dependence on these parameters.
      The dataset 1a was used to evaluate the dependency of the Mu.f, Sigma.f and WAVE.f values as a function of slice thickness, kernel and reconstruction algorithm.
      For comparison, this variability was assessed on another lung metric, previously described by several authors [

      M. Mascalchi G. Camiciottoli S. Diciotti Lung densitometry: Why, how and when J Thorac Dis 9 2017 3319 45 https://doi.org/10.21037/jtd.2017.08.17.

      ,

      Lynch D. Progress in Imaging COPD, 2004 - 2014. Chronic Obstr Pulm Dis J COPD Found 2014. https://doi.org/10.15326/jcopdf.1.1.2014.0125.

      ,
      • Gattinoni L.
      • Pesenti A.
      • Bombino M.
      • Baglioni S.
      • Rivolta M.
      • Rossi F.
      • et al.
      Relationships Between Lung Computed Tomographic Density, Gas Exchange, and PEEP in Acute Respiratory Failure.
      ,
      • Chen A.
      • Karwoski R.A.
      • Gierada D.S.
      • Bartholmai B.J.
      • Koo C.W.
      Quantitative CT analysis of diffuse lung disease.
      ,

      Lanza E, Muglia R, Bolengo I, Santonocito OG. Quantitative Chest CT analysis in COVID-19 to predict the need for oxygenation support and intubation. PrePrint 2020:1–18. https://doi.org/10.21203/rs.3.rs-30481/v1.

      ], calculated integrating histogram within fixed threshold in the HU range −950 ÷ -700, hereafter called WAVE.th. Dataset 2 allowed to investigate the correlation between the Gaussian curve metrics and the respiratory phase. Different level of inspiration can affect the central value (Mu.f) of the Gaussian function but the area under this curve should represent an estimation of the well-aerated lung volume, independently of the respiratory phase.
      The tidal volume (TD) was calculated as the maximum difference in lung volume during the respiratory cycle and patients with TD lower than 390 mL were excluded from dataset 2 [
      • Ekström M.
      • Sundh J.
      • Schiöler L.
      • Lindberg E.
      • Rosengren A.
      • Bergström G.
      • et al.
      Absolute lung size and the sex difference in breathlessness in the general population.
      ,
      • Needham C.D.
      • Rogan M.C.
      • McDonald I.
      Normal Standards for Lung Volumes, Intrapulmonary Gas-mixing, and Maximum Breathing Capacity.
      ]. With respect to tumor location, contralateral lungs of the remaining patients, were analyzed and the metrics obtained at the various respiratory phases were correlated with the percentage variation of the lung volume.

      2.3.3 Biomarkers in COVID-19 patients imaging

      For dataset 3, in addition to all the metrics described above (Mu.f, Sigma.f, WAVE.f, WAVE.th), we calculated also standard first-order radiomic features (i.e. average HU, Skewness, Kurtosis and 0.25, 0.5, 0.75 and 0.90 percentiles) from the CT histogram of the whole lung volume. For comparison of healthy and affected lungs, the same biomarkers were calculated also in patients of dataset 1.
      A new biomarker was introduced to assess local change in lung density with respect to the healthy parenchyma and defined as:
      ΔHUMu-Avgi=HU¯i-Mu.f
      (3)


      In Eq. (3), Mu.f is the fit parameter relative of the histogram for the entire lung whereas HU¯i is the average value of voxels in a specific region i of the lung. This biomarker was then calculated in each of the 24 lung subregions of COVID-19 patients.
      For each patient, the range and the median values over all sub regions were considered as a metric related to the severity of the disease. We calculated the average and the standard deviation values over the 20 patients.

      2.4 Data analysis

      Statistical analysis was performed by using the Real Statistics Resource Pack software (release 6.8, www.realstatistics.com). Saphiro-Wilk and Levene tests were used to assess the normality of the distributions and the homogeneity of the variance.

      2.4.1 Dataset 1

      The parameters of the Gaussian function were calculated with both the polynomial and the exponential equations, fitting data points of histogram calculated with 5 HU using a 30% and 90% of the modal value as starting and ending points. The two models were compared using the relationship of Eq. (2). Differences between Height.f, Mu.f, Sigma.f and WAVE.f values resulting from the polynomial and exponential models were assessed for each healthy subject of dataset 1 and their correlation was assessed using the Pearson coefficient.
      To estimate the statistical significance of the differences in WAVE.f calculation due to the bin width and to the range extremes for the Gaussian fitting, analysis of variance (ANOVA) was performed (significance level of 0.05), changing one parameter at a time and keeping the other two fixed.

      2.4.2 Dataset 1a

      ANOVA was applied (significance level of 0.05) to the Mu.f, Sigma.f, WAVE.f and WAVE.th values calculated in CT images reconstructed from the same sinogram changing the reconstruction parameters (slice thickness, reconstruction algorithm, reconstruction kernel).
      Paired t-tests were then calculated between metrics calculated from the VCT (3 mm as slice thickness, IR-S1 and Bl57 as reconstruction algorithm and kernel) and the other post-processed series.
      A p value<0.05 was considered indicating a significant difference.

      2.4.3 Dataset 2

      Pearson correlation coefficient was calculated to estimate the linear correlation between the metrics described (Mu.f, Sigma.f, WAVE.f, WAVE.th) and the lung volume, expressed as a percentage of the maximum volume, at the 10 different phases of the respiratory cycle. According to the distribution of the Pearson correlation coefficient with n = 10, values above 0.632 and 0.765 were considered indicating significant and strongly significant correlation, respectively (significance level: 0.05 and 0.01). McNemar test was used to assess the significance of the differences in the number of patients with WAVE.f and WAVE.th correlated with the respiratory phase.

      2.4.4 Dataset 3

      A t-test was used to compare results of this dataset with those obtained in dataset 1 (significance level of 0.05). A description of lung disease is given by a graphical representation of the CT lung histogram metrics in the 24 lung subregions, hereafter called lungogram.

      3. Results

      The properties and results of the metrics and biomarkers (WAVE.f, WAVE.th and ΔHUMu-Avg) were studied in relation to the specific datasets described above.

      3.1 Dataset 1

      Average differences for Height.f, Mu.f, Sigma.f and WAVE.f calculated with polynomial and exponential model were 0.1%, 0.2 HU, 0.3 HU and 0.3%, respectively. Pearson coefficients for these metrics were 1.00, 1.00, 0.979 and 0.988. Since CT histogram are generally reported in linear scale, hereafter, only results of the exponential fitting will be reported.
      Fitting parameters, such as bin width, starting point and end point of the fitting range, did not affected significantly the WAVE.f calculation. According to the ANOVA tests, the p-value was>0.812 for all them.
      In our analysis, a value of 5 HU as bin width was selected since it provided a trade-off between smoothness and resolution in the histogram. In healthy subjects, the choice of starting and ending point is not relevant, since CT histogram of the healthy parenchyma does not deviate importantly from a Gaussian distribution. In this work, considering the application on compromised lungs, 30% on the left side of the modal histogram value was set as starting point (see the point (a) of Fig. 2). In this way, the fitting range should avoid the less dense part of the histogram possibly due to the presence of emphysema (Fig. 2). As ending point, 90% on the right side of the modal histogram value was set (see the point (b) of Fig. 2). This choice should reduce the use of the points on the right side of the histogram in patients with increased lung density, as occurs in particular for the COVID-19 dataset.

      3.1.1 Dataset 1a

      ANOVA found statistically significant differences changing slice thickness and iterative strength for the estimation of WAVE.f (p < 0.001, p < 0.001), WAVE.th (p < 0.001, p < 0.001) and Sigma.f (p = 0.026, p < 0.001). For Mu.f the differences were not statistically significant (p = 0.629, p = 0.896). All results regarding analysis of dataset 1a are summarized in Table 3.
      Table 3Mean and Standard Deviation of Mu.f, Sigma.f, WAVE.f and WAVE.th values calculated for the dataset 1a, changing the reconstruction parameters (slice thickness, reconstruction algorithm, kernel).
      Kernel and reconstruction algorithmSlice ThicknessMu.fSigma.fWAVE.fWAVE.th
      (mm)(HU)(HU)(%)(%)
      Bl57 IR-11 mm−887(±26)99(±10)82(±2)66(±3)
      Bl57 IR-13 mm−868(±26)74(±11)84(±2)75(±4)
      Bl57 IR-15 mm−867(±25)62(±7)80(±2)77(±3)
      Bl57 FBP3 mm−870(±26)83(±11)83(±2)72(±3)
      Bl57 IR-13 mm−868(±26)74(±11)84(±2)75(±3)
      Bl57 IR-33 mm−866(±25)51(±7)80(±3)82(±3)
      Bl57 IR-53 mm−864(±24)34(±6)77(±4)87(±2)
      Bl57 IR-13 mm−868(±26)74(±11)84(±2)75(±4)
      Br38 IR-13 mm−852(±25)30(±7)72(±6)88(±3)
      The Gaussian fit parameter Mu.f did not showed significant differences in paired t-test between the standard reconstruction (VCT series) and the other series comparing 3 and 5 mm as slice thickness or IR-S1 and IR-S3 as reconstruction algorithm (Table 4). Sigma.f and WAVE.th had always significant differences changing the reconstruction parameters, while WAVE.f did not differ using slice thickness of 1 or 3 mm.
      Table 4P-values results of paired t-test on dataset 1a changing the reconstruction parameters (slice thickness, reconstruction algorithm, kernel), keeping as reference the standard series (3 mm, IR-S1, Bl57).
      Mu.fSigma.fWAVE.fWAVE.th
      Slice Thickness1 vs 3 mm<0.001<0.0010.09<0.001
      1 vs 5 mm<0.001<0.001<0.001<0.001
      3 vs 5 mm0.8<0.001<0.001<0.001
      Reconstruction AlgorithmFBP vs IR-S10.04<0.0010.01<0.001
      FBP vs IR-S3<0.001<0.001<0.001<0.001
      IR-S1 vs IR-S30.2<0.001<0.001<0.001
      IR-S1 vs IR-S50.03<0.001<0.001<0.001
      KernelBl57 vs Br38<0.001<0.001<0.001<0.001

      3.2 Dataset 2

      From the second cohort, 5 patients were excluded for their limited difference in lung volume in the respiratory cycle (TD < 390 mL). For the remaining 15 patients, the results of the correlation between the calculated metrics and lung expansion are reported in Table 5. WAVE.f resulted correlated significantly in 3 and strongly significant in 1 out of 15 cases. As expected, significant correlation was found for Mu.f in all cases. Sigma.f resulted significantly correlated in 11 cases. For WAVE.th, the correlation was strongly significant for 12 patients.
      Table 5Number of cases with lung metrics calculated in 4DCT images, correlated with respiratory phase.
      AlphaPearson coefficientMu.fSigma.fWAVE.fWAVE.th
      0.05 (significant correlation)> 0.632n1511312
      0.01 (strongly significant correlation)> 0.765n148112
      McNemar test applied to the number of cases reported in Table 5 for WAVE.f and WAVE.th returned p-values of 0.009 and 0.001 for significant and strongly significant correlation with respiratory cycle indicating a significant difference between the two metrics in the correlation with the respiratory cycle.
      The effect of respiratory cycle on the lung metrics calculated in CT images for two patients (P104, P111) from the dataset 2 representative of results of Table 5 is reported in Fig. 3. Linear regression lines are displayed when a strongly significant correlation was found between lung volume and lung biomarkers (Mu.f, Sigma.f and WAVE.f and WAVE.th).
      Figure thumbnail gr3
      Fig. 3The effect of respiratory cycle on the lung metrics calculated in CT images for two patients from the dataset 2 representative of results reported in . (A) and (B): CT lung relative histogram of the two patients in minimum and maximum inhale phases. (C) and (D): fit parameters Mu.f and Sigma.f are plotted versus the lung volume. (E) : lung volume plotted versus the respiratory phase. (F): WAVE metrics plotted versus lung volume. Linear regression lines are displayed when strong correlation was found.

      3.3 Dataset 3

      Metrics calculated for the whole lung in COVID-19 patients of dataset 3 were compared with results of normal lung patients of dataset 1. Boxplot of the lung metrics calculated in dataset 1 (healthy lung) and dataset 3 (COVID-19 patients) are reported in Fig. 4. Differences were significant for all the metrics. A p-value of 0.007 was found for Mu.f comparison while for all the other metrics p-values were < 0.001. In healthy subjects average (standard deviation) values for WAVE.f and WAVE.th resulted 84% (2%) and 75% (4%). Results of the biomarker ΔHUMu-Avg calculated in the 24 subregions for the dataset 3 are summarized in Table 6. The average (standard deviation) values for range and median of ΔHUMu-Avg were, respectively, 411 (127) HU and 179 (59) HU. Examples of analysis performed for two patients, including histogram, Gaussian fit and lungogram, are reported in Fig. 5. The difference between CT lung histogram and Gaussian fit data, indicated with the blue line, represents the distribution of all non-healthy lung parenchyma like vessels, airways and disease-involved tissue.
      Figure thumbnail gr4
      Fig. 4Boxplot of the lung metrics calculated in dataset 1 (normal lung) and dataset 3 (COVID-19 patients).
      Table 6Minimum, maximum, median values and range of the biomarker ΔHUMu-Avg calculated over the 24 subregions for all patients of the dataset 3. The last two rows report the average and the standard deviation over the all dataset.
      Patient IDΔHUMu-Avg(HU)
      MinimumMaximumRangeMedian
      P_01−323123476
      P_0287301214165
      P_0355528472299
      P_0431422390202
      P_0526293267120
      P_0628331304155
      P_0741666625260
      P_0818352334115
      P_0973372299144
      P_10−11467479155
      P_1160617557215
      P_1226339313144
      P_1314647633117
      P_1414472458177
      P_1549386337216
      P_1642480438175
      P_1730522493197
      P_1854644590284
      P_1949396347232
      P_2014446432144
      Average35446411179
      Standard Deviation2512712759
      Figure thumbnail gr5
      Fig. 5Examples of patients analyzed with the methods described in the materials and methods section. Segmented lungs are displayed in axial and coronal views. Local changes in density are reported in lungograms. In the plots, black curve represents the CT relative histogram of the entire lung. Yellow dots are the data points selected to calculate Gaussian fit (red curve). The difference between CT lung relative histogram and Gaussian fit is reported in blue while the HU range used for WAVE.th calculations is reported by two green vertical lines. The relative metrics are reported with the same colors as the curves in the graphs on which they are calculated. Lung histogram features (HU average, skewness, kurtosis, and percentiles) are reported in black, Gaussian fit parameters (Mu.f, Sigma.f) and WAVE.f are reported in red while WAVE.th is reported in green. Case A: example of healthy lung from dataset 1. The Gaussian fit is well overlapped to the CT relative histogram, with a slight deviation on the right side of the black curve, caused likely by the inclusion in the lung segmentation of major vessels. Case B: example of patient with cancer in the left lung, from dataset 2. The presence of the tumor mass and the consequent distortion of the surrounding structures increase the deviation between the Gaussian fit and the CT relative histogram. Case C: COVID-19 patient with an extended ground glass opacification of the lung that enlarges the Gaussian fit. In addition to the fact that Mu.f value is on the right side of the green window (-950/-700 HU), it causes an important difference between the two WAVE metrics. The lungogram shows that the higher values of the biomarker ΔHUMu-Avg are found in the dorsal regions. Case D: Widespread disease with only some portions of healthy lung, but sufficient to produce a suitable histogram for Gaussian fit. In this case, WAVE.f, and WAVE.th calculations returned similar results even if Mu.f value is not in the central position of the range used to calculate WAVE.th. The lungogram shows that the higher values of ΔHUMu-Avg are found in the dorsal regions. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

      4. Discussion

      In this work a patient-specific and automatic model, based on a Gaussian fit of relative CT lung histogram, was described to quantify aerated and pathologic lung regions in healthy controls and in patients with interstitial pneumonia caused by SARS-COV-2. The parameters Mu.f and Sigma.f, which characterize the Gaussian model, are related to lung inflation and noise in the images of the healthy parenchyma, resepctively. WAVE.f is a new metric suitable for quantification of the well-aerated lungs that showed lower variability with physiological and technical factors than others already proposed.
      CT lung histogram fitting has already been reported in literature and metrics derived from histogram analysis have been used in the past. Obert et Al.[
      • Obert M.
      • Kampschulte M.
      • Limburg R.
      • Barańczuk S.
      • Krombach G.A.
      Quantitative computed tomography applied to interstitial lung diseases.
      ] proposed a model in which cumulative CT lung histogram data were fitted with a logistic growth function in order to classify normal and pathological lungs. Despite the good results, this logistic model had a principle mathematical approach. In contrast, the model we propose takes into account the physical properties of the pulmonary parenchyma and the nature of noise of CT images. In this work, we used a non-linear equation to calculate the Mu.f and Sigma.f since Gaussian fit is readily available within the program used to develop our model. Alternatively, a second order polynomial fitting of the logarithm of the histogram frequencies could be used making the procedure more independent from the platform used. However, the difference between the WAVE.f results was negligible for clinical purposes and the exponential fitting offers the advantage of having a clear graphical interpretation of CT lung histogram.
      Usually, quantitative lung analysis applies fixed HU thresholds in order to divide lung volume in compartments with different aerations. However, a fixed HU range does not take in account the inter-patient variability or the phase of the respiratory cycle. This fact is highlighted in the examples reported in Fig. 5, where the Gaussian peak that accurately fits data points around the histogram modal value is close to the opposite limits of the fixed range (Mu.f = −765 HU for Case C, Mu.f =−911 HU for Case D). The inter-patient variability of the histogram modal value, which results in this asymmetry with respect to the fixed HU range, can be ascribed to clinical or physiological factors which only a patient-specific metric can take into account.
      The model of this study was tested in two cohorts of healthy and critically ill patients, representing two opposite conditions of lung status. Under the hypothesis formulated here, a Gaussian curve describes the radiological aspects of healthy lung tissue: the central value represents the average density of the well-aerated tissue and the width includes all sources of image noise due to reconstruction parameters and anatomical texture. The two model parameters Mu.f and Sigma.f inevitably depend on the histogram data points used for fitting but the resulting WAVE.f values are independent of the fitting range used, at least in healthy lungs. Even in cases with severe lung impairment due to interstitial pneumonia, the modal value of the histograms was representative of healthy lung tissue and the Gaussian fit properly represents the estimation of the ventilated lung as it takes into account the intra-patient variability. However, in a limited fraction of subjects having most of the lung volume compromised by opacities or solidifications, the histogram modal value of the whole lung can results with unusual values, not included in the range typically characterizing the health lung tissue, i.e. between −950 and −700 HU. In these cases, the Gaussian model may not be directly applicable as the fit on the selected histogram data could return for Mu.f values that are not representative of the healthy tissue.
      To take this into account, the definition of the fit range should not consider the modal value of the whole lung histogram but the presence of relative maxima in a limited HU range. This is evident in the example shown in Fig. 6, where the modal value of the histogram is greater than −700HU. The value of Mu.f obtained fitting data points around the modal value, or even in a limited range as shown in Fig. 6-A, is higher than −700 underlining the non-applicability of the model under these conditions.
      Figure thumbnail gr6
      Fig. 6Example of a COVID-19 patient with widespread opacification, displayed in axial and sagittal views. Plot A: the modal value in the CT lung histogram is greater than −700 HU and the definition of the fit range based on this value is not correct for the Gaussian model. Plot B: Data points for Gaussian fit are defined by the local maximum identified by first derivative. Plot C: Analysis performed in the anterior-medial of the central portion of the right lung where, as showed by the yellow cross-hair in the CT images, normal lung is prevalent. The light red shadowed areas under the curves represent the lung volume classified as “well aerated” according to the Gaussian model but not according to theWAVE.th definition. Analogously, the light blue shadowed areas under the curves represent the lung volume classified as “well aerated” according to the fixed threshold model but not according to the Gaussian metric. The difference between the two metrics arise when these two areas do not offset each other. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
      On the other hand, in the case of patients with emphysema [
      • Schroeder J.D.
      • McKenzie A.S.
      • Zach J.A.
      • Wilson C.G.
      • Curran-Everett D.
      • Stinson D.S.
      • et al.
      Relationships Between Airflow Obstruction and Quantitative CT Measurements of Emphysema, Air Trapping, and Airways in Subjects With and Without Chronic Obstructive Pulmonary Disease.
      ], because of the phenomenon known as “air trapping”, the Mu.f values might be underestimated if the fit range values are not adequate. To overcome the limitations of these specific cases, it is possible to select the histogram points to be fitted in a more tailored way, for example using additional conditions on the first or second derivative (Fig. 6-B). Alternatively, useful information could come from histograms of sub-regions where the disease is less prevalent and the modal value is representative of the healthy tissue, such as in the sub-regions of the same patient in Fig. 6-C where the Gaussian model can be applied.
      Regarding the robustness, the metrics for QCT assessment should be independent from technical and physiological bias: analysis of dataset 1a and dataset 2 was aimed at studying the relationship of the metrics with respiratory cycle and with image reconstruction parameters.
      The sources of variability in QCT for lung analysis are well known and reported in literature [

      Mackin D, Fave X, Zhang L, Fried D, Yang J, Brian Taylor, et al. Measuring computed tomography scanner variability of radiomics features. Invest Radiol 2015;50:757–65. https://doi.org/10.1097/RLI.0000000000000180.

      ,
      • Reske A.W.
      • Busse H.
      • Amato M.B.P.
      • Jaekel M.
      • Kahn T.
      • Schwarzkopf P.
      • et al.
      Image reconstruction affects computer tomographic assessment of lung hyperinflation.
      ,
      • Chen-Mayer H.H.
      • Fuld M.K.
      • Hoppel B.
      • Judy P.F.
      • Sieren J.P.
      • Guo J.
      • et al.
      Standardizing CT lung density measure across scanner manufacturers.
      ]. A protocol standardization is recommended to reduce results variability, but it is achievable only in prospective studies. Neither WAVE.f nor WAVE.th metrics resulted completely independent from all the reconstruction parameters. However, no significant differences were found in WAVE.f values calculated in images with slice thickness in the range 1–3 mm, the most used values for high resolution lung CT protocols.
      An increase of slice thickness, as well as an increase of the strength of iterative algorithm, implies a decrease in image noise magnitude that impacts systematically on the CT lung histograms and on the derived metrics. Comparing WAVE.f and WAVE.th values calculated in images reconstructed with different parameters, an overall lower variability was observed in the Gaussian model’s metric due to its intrinsic ability to fit the actual data in the histogram.
      A significant difference was found in all comparison when different algorithms were selected. The nonlinear effects of SAFIRE algorithm may explain the increase of differences found with the increase of the IR strength [

      Solomon J, Samei E. Quantum noise properties of CT images with anatomical textured backgrounds across reconstruction algorithms: FBP and SAFIRE. Med Phys 2014;41. https://doi.org/10.1118/1.4893497.

      ]. Nevertheless, the most limited difference in the WAVE.f results was found comparing a moderate strength of iteration and FBP algorithm. For this reason, VCT series, using Bl57 IR-S1 and 3 mm as slice thickness, was considered suitable also for quantitative purpose.
      Replacing the Bl57 kernel with the medium smooth Br38 kernel, not generally used for diagnostic purposes in chest imaging, implies the greatest differences in WAVE.f. However, we added this reconstruction series since images in public dataset 2 were available only with a soft kernel but medium-smooth kernel are not generally used for VCT of the lung parenchyma.
      Another challenging task in medical imaging is organ motion [
      • Nehmeh S.A.
      • Erdi Y.E.
      • Ling C.C.
      • Rosenzweig K.E.
      • Squire O.D.
      • Braban L.E.
      • et al.
      Effect of respiratory gating on reducing lung motion artifacts in PET imaging of lung cancer.
      ,
      • Liu C.
      • Pierce L.A.
      • Alessio A.M.
      • Kinahan P.E.
      The impact of respiratory motion on tumor quantification and delineation in static PET/CT imaging.
      ,
      • Rietzel E.
      • Chen G.T.Y.
      • Choi N.C.
      • Willet C.G.
      Four-dimensional image-based treatment planning: Target volume segmentation and dose calculation in the presence of respiratory motion.
      ]. In particular, for quantitative imaging assessment, the differences in lung volume due to the respiratory cycle result in different values of tissue density. This is clearly visible in Fig. 3 and in Table 5, where the correlation for Mu.f was always significant. Sigma.f also showed a significant correlation in most of the cases, affecting the values of WAVE.th that, as consequence, showed a strongly significant correlation with lung inflation in 12 out of 15 of analyzed cases. On the other hand, WAVE.f resulted more stable with respiratory phases than other metrics. This is due to the customized properties of the proposed model that takes into account the inter- and intra- subject variability. It must be stressed that these results refer to radiotherapy patients, specifically trained to follow a shallow and regular breathing, while during diagnostic chest examination, the CT scan is performed by asking the patient to have deep breath, compatibly with his pathological state.
      The estimation of the fractional volume of well aerated lung is calculated for both WAVE metrics from the integration of the CT lung histogram but with a substantial methodological difference: while the WAVE.th metric uses fixed HU range to identify the histogram area corresponding to the healthy lung parenchyma, in the proposed method the integration range is chosen according to an image-specific model. These different approaches are represented in Fig. 6-B and 6-C. The light red shadowed areas under the curves represent the lung volume classified as “well aerated” according to the Gaussian model but not according to the WAVE.th definition. Analogously, the light blue shadowed areas under the curves represent the lung volume classified as “well aerated” according to the fixed threshold model but not according to the Gaussian metric. When the differences between these volumes are not compensated, a difference between the two WAVE metrics occurs. The relationship between WAVE.f and WAVE.th is shown in Fig. 7. As expected, lung opacifications and solidifications in COVID-19 patients reduce WAVE values for both metrics, with WAVE.th systematically lower than WAVE.f. Moreover, in healthy subjects no correlation exists because the variability of WAVE.th, that do not consider the actual respiratory phase, is not hidden by the inter-patient disease variability.
      Figure thumbnail gr7
      Fig. 7Scatterplot of WAVE.f and WAVE.th values for patients in dataset 1 (healthy lung) and dataset 3 (COVID-19 patients).
      In healthy subjects, an average value of 84% (range: 79%-86%) was found for WAVE.f. This result represents the percentage of healthy parenchymal tissue of the entire lung volume in non-pathological lung images. Although the assessment of the accuracy of WAVE metrics was not amongst the aims of this work, a correspondence was found with results of morphometric studies. Townsley reported an average value of 84% (range: 77–87%) as fraction of overall anatomic lung defined as parenchyma [
      • Townsley M.I.
      Structure and composition of pulmonary arteries, capillaries, and veins.
      ]. By contrast, WAVE.th calculated in the same cohort of healthy subjects, showed lower results and a higher variability.
      In the comparison between results of datasets 1 and 3, for all the metrics the t-test showed significant differences. The discrepancy between the two cohorts found for Mu.f can be explained as the limited capabilities of deep breath holding in patients with severe lung impairments that result in higher density values of parenchymal tissue.
      The other radiomic features calculated from histograms (HU mean, Skewness and Kurtosis) have similar trends as reported in literature to distinguish healthy and pathological lungs [
      • Ash S.Y.
      • Harmouche R.
      • Vallejo D.L.L.
      • Villalba J.A.
      • Ostridge K.
      • Gunville R.
      • et al.
      Densitometric and local histogram based analysis of computed tomography images in patients with idiopathic pulmonary fibrosis.
      ,
      • Best A.C.
      • Meng J.
      • Lynch A.M.
      • Bozic C.M.
      • Miller D.
      • Grunwald G.K.
      • et al.
      Idiopathic pulmonary fibrosis: Physiologic tests, quantitative CT indexes, and CT visual scores as predictors of mortality.
      ].
      As expected, also WAVE metrics clearly discerned the two cohorts, but its values, unlike the first-order radiomic features, can give quantitative information about the well-aerated lung that can be useful in the management of patients with ARDS [
      • Gattinoni L.
      • Caironi P.
      • Cressoni M.
      • Chiumello D.
      • Ranieri V.M.
      • Quintel M.
      • et al.
      Lung recruitment in patients with the acute respiratory distress syndrome.
      ].
      Moreover, WAVE.f values are referred to the entire lung but interstitial pneumonia from SARS-CoV-2 produces solidifications and opacifications distributed heterogeneously in the lung parenchyma with different level of severity. Therefore, local changes in lung density can be assessed calculating ΔHUMu-Avg in any specific lung region. A graphical representation of this biomarker in the lungogram provides spatial and quantitative information of the patient's lung status in a single 2D picture of easy interpretation and suitable for clinical decision or inter patients’ comparison.
      This work has some limitations. First of all, the model is not suitable for all cases. When the disease affects most of the lung volume and the modal value of the CT numbers is shifted towards unexpected high values, a different approach must be used to extend the applicability of the model to the entire lung. Possible solutions include different criteria for the definition of the fit range or taking into account the histogram of the lung sub-regions. Furthermore, the evaluation of the impact of the reconstruction parameters on lung metrics was performed using only one scanner. However, such a study regarding both reconstruction and acquisition parameters is feasible only with physical phantom or with a very high number of cases involved.

      5. Conclusion

      A method to analyse CT lung images based on the Gaussian fit of the histogram data has been developed and characterized.
      In healthy lungs, WAVE.f, a new quantitative metric derived from physics assumptions and with physiological significance, demonstrates lower dependencies from technical or physiological parameters with respect to the already reported equivalent metrics and its values were in good agreement with morphometric studies.
      The complex patterns of lung diseases, such as those resulting from SARS-CoV-2 pneumonia, can be described by appropriate metrics calculated locally. The biomarker ΔHUMu-Avg is an absolute measure in Hounsfield Units of lung density and its values calculated in 24 lung subregions of COVID-19 patients combines quantitative and spatial information.
      Finally, a validation of WAVE metrics is mandatory before its use for clinical decision. A future work using a larger sample of clinical images and functional data can be addressed to verify the hypothesis on which this model is built and to assess accuracy of the WAVE.f.

      References

        • Zhou F.
        • Yu T.
        • Du R.
        • Fan G.
        • Liu Y.
        • Liu Z.
        • et al.
        Clinical course and risk factors for mortality of adult inpatients with COVID-19 in Wuhan, China: a retrospective cohort study.
        Lancet. 2020; 395: 1054-1062https://doi.org/10.1016/S0140-6736(20)30566-3
        • Chung M.
        • Bernheim A.
        • Mei X.
        • Zhang N.
        • Huang M.
        • Zeng X.
        • et al.
        CT imaging features of 2019 novel coronavirus (2019-NCoV).
        Radiology. 2020; 295: 202-207https://doi.org/10.1148/radiol.2020200230
        • Rubin G.D.
        • Ryerson C.J.
        • Haramati L.B.
        • Sverzellati N.
        • Kanne J.P.
        • Raoof S.
        • et al.
        The Role of Chest Imaging in Patient Management During the COVID-19 Pandemic.
        Chest. 2020; 158: 106-116https://doi.org/10.1016/j.chest.2020.04.003
        • Ghetti C.
        • Ortenzia O.
        • Maddalo M.
        • Altabella L.
        • Sverzellati N.
        Dosimetric and radiation cancer risk evaluation of high resolution thorax CT during COVID-19 outbreak.
        Phys Medica. 2020; 80: 119-124https://doi.org/10.1016/j.ejmp.2020.10.018
        • Salehi S.
        • Abedi A.
        • Balakrishnan S.
        • Gholamrezanezhad A.
        Coronavirus disease 2019 (COVID-19) imaging reporting and data system (COVID-RADS) and common lexicon: a proposal based on the imaging data of 37 studies.
        Eur Radiol. 2019; 2020: 2019https://doi.org/10.1007/s00330-020-06863-0
        • Hansell D.M.
        • Bankier A.A.
        • MacMahon H.
        • McLoud T.C.
        • Müller N.L.
        • Remy J.
        Fleischner Society: Glossary of terms for thoracic imaging.
        Radiology. 2008; 246: 697-722https://doi.org/10.1148/radiol.2462070712
      1. S. Simpson F.U. Kay S. Abbara S. Bhalla J.H. Chung M. Chung et al. Radiological Society of North America Expert Consensus Statement on Reporting Chest CT Findings Related to COVID-19. Endorsed by the Society of Thoracic Radiology, the American College of Radiology, and RSNA. Radiol Cardiothorac Imaging 2020;2:e200152. 10.1148/ryct.2020200152.

        • Li K.
        • Fang Y.
        • Li W.
        • Pan C.
        • Qin P.
        • Zhong Y.
        • et al.
        CT image visual quantitative evaluation and clinical classification of coronavirus disease (COVID-19).
        Eur Radiol. 2020; 30: 4407-4416https://doi.org/10.1007/s00330-020-06817-6
        • Yang R.
        • Li X.
        • Liu H.
        • Zhen Y.
        • Zhang X.
        • Xiong Q.
        • et al.
        Chest CT Severity Score: An Imaging Tool for Assessing Severe COVID-19.
        Radiol Cardiothorac Imaging. 2020; 2https://doi.org/10.1148/ryct.2020200047
        • Sumikawa H.
        • Johkoh T.
        • Yamamoto S.
        • Yanagawa M.
        • Inoue A.
        • Honda O.
        • et al.
        Computed tomography values calculation and volume histogram analysis for various computed tomographic patterns of diffuse lung diseases.
        J Comput Assist Tomogr. 2009; 33: 731-738https://doi.org/10.1097/RCT.0b013e31818da65c
      2. M. Mascalchi G. Camiciottoli S. Diciotti Lung densitometry: Why, how and when J Thorac Dis 9 2017 3319 45 https://doi.org/10.21037/jtd.2017.08.17.

        • Wang J.
        • Li F.
        • Li Q.
        Automated segmentation of lungs with severe interstitial lung disease in CT.
        Med Phys. 2009; 36: 4592-4599https://doi.org/10.1118/1.3222872
        • Park B.
        • Park H.
        • Lee S.M.
        • Seo J.B.
        • Kim N.
        Lung Segmentation on HRCT and Volumetric CT for Diffuse Interstitial Lung Disease Using Deep Convolutional Neural Networks.
        J Digit Imaging. 2019; 32: 1019-1026https://doi.org/10.1007/s10278-019-00254-8
        • Nakagawa H.
        • Nagatani Y.
        • Takahashi M.
        • Ogawa E.
        • Van Tho N.
        • Ryujin Y.
        • et al.
        Quantitative CT analysis of honeycombing area in idiopathic pulmonary fibrosis: Correlations with pulmonary function tests.
        Eur J Radiol. 2016; 85: 125-130https://doi.org/10.1016/j.ejrad.2015.11.011
      3. Lynch D. Progress in Imaging COPD, 2004 - 2014. Chronic Obstr Pulm Dis J COPD Found 2014. https://doi.org/10.15326/jcopdf.1.1.2014.0125.

        • Ash S.Y.
        • Harmouche R.
        • Vallejo D.L.L.
        • Villalba J.A.
        • Ostridge K.
        • Gunville R.
        • et al.
        Densitometric and local histogram based analysis of computed tomography images in patients with idiopathic pulmonary fibrosis.
        Respir Res. 2017; 18: 1-11https://doi.org/10.1186/s12931-017-0527-8
        • Gattinoni L.
        • Caironi P.
        • Pelosi P.
        • Goodman L.R.
        What Has Computed Tomography Taught Us about the Acute Respiratory Distress Syndrome?.
        Am J Respir Crit Care Med. 2001; 164: 1701-1711https://doi.org/10.1164/ajrccm.164.9.2103121
        • Gattinoni L.
        • Pesenti A.
        • Bombino M.
        • Baglioni S.
        • Rivolta M.
        • Rossi F.
        • et al.
        Relationships Between Lung Computed Tomographic Density, Gas Exchange, and PEEP in Acute Respiratory Failure.
        Anesthesiology. 1988; 69: 824-832https://doi.org/10.1097/00000542-198812000-00005
        • Nishiyama A.
        • Kawata N.
        • Yokota H.
        • Sugiura T.
        • Matsumura Y.
        • Higashide T.
        • et al.
        A predictive factor for patients with acute respiratory distress syndrome: CT lung volumetry of the well-aerated region as an automated method.
        Eur J Radiol. 2020; 122https://doi.org/10.1016/j.ejrad.2019.108748
        • Gattinoni L.
        • Caironi P.
        • Cressoni M.
        • Chiumello D.
        • Ranieri V.M.
        • Quintel M.
        • et al.
        Lung recruitment in patients with the acute respiratory distress syndrome.
        N Engl J Med. 2006; https://doi.org/10.1056/NEJMoa052052
        • Chen A.
        • Karwoski R.A.
        • Gierada D.S.
        • Bartholmai B.J.
        • Koo C.W.
        Quantitative CT analysis of diffuse lung disease.
        Radiographics. 2020; 40: 28-43https://doi.org/10.1148/rg.2020190099
      4. D. Colombi F.C. Bodini M. Petrini G. Maffi N. Morelli G. Milanese et al. Well-aerated Lung on Admitting Chest CT to Predict Adverse Outcome in COVID-19 Pneumonia Radiology 2020:201433. 10.1148/radiol.2020201433.

        • Matsuoka S.
        • Yamashiro T.
        • Matsushita S.
        • Kotoku A.
        • Fujikawa A.
        • Yagihashi K.
        • et al.
        Quantitative CT evaluation in patients with combined pulmonary fibrosis and emphysema: Correlation with pulmonary function.
        Acad Radiol. 2015; 22: 626-631https://doi.org/10.1016/j.acra.2015.01.008
        • Ohkubo H.
        • Kanemitsu Y.
        • Uemura T.
        • Takakuwa O.
        • Takemura M.
        • Maeno K.
        • et al.
        Normal lung quantification in usual interstitial pneumonia pattern: The impact of threshold-based volumetric CT analysis for the staging of idiopathic pulmonary fibrosis.
        PLoS ONE. 2016; 11: 1-13https://doi.org/10.1371/journal.pone.0152505
      5. Lanza E, Muglia R, Bolengo I, Santonocito OG. Quantitative Chest CT analysis in COVID-19 to predict the need for oxygenation support and intubation. PrePrint 2020:1–18. https://doi.org/10.21203/rs.3.rs-30481/v1.

        • Chiumello D.
        • Langer T.
        • Vecchi V.
        • Luoni S.
        • Colombo A.
        • Brioni M.
        • et al.
        Low-dose chest computed tomography for quantitative and visual anatomical analysis in patients with acute respiratory distress syndrome.
        Intensive Care Med. 2014; 40: 691-699https://doi.org/10.1007/s00134-014-3264-1
        • Cressoni M.
        • Gallazzi E.
        • Chiurazzi C.
        • Marino A.
        • Brioni M.
        • Menga F.
        • et al.
        Limits of normality of quantitative thoracic CT analysis.
        Crit Care. 2013; 17: R93https://doi.org/10.1186/cc12738
        • Traverso A.
        • Wee L.
        • Dekker A.
        • Gillies R.
        Repeatability and reproducibility of radiomic features: a systematic review.
        Int J Radiat Oncol Biol Phys. 2018; 102: 1143-1158https://doi.org/10.1016/j.ijrobp.2018.05.053
        • Newell Jr., J.D.
        • Tschirren J.
        • Peterson S.
        • Beinlich M.
        • Sieren J.
        Quantitative CT of interstitial lung disease.
        Semin Roentgenol. 2019; 54: 73-79https://doi.org/10.1053/j.ro.2018.12.007
        • Tian Y.
        • Miao J.
        • Liu Z.
        • Huang P.
        • Wang W.
        • Wang X.
        • et al.
        Availability of a simplified lung ventilation imaging algorithm based on four-dimensional computed tomography.
        Phys Medica. 2019; 65: 53-58https://doi.org/10.1016/j.ejmp.2019.08.006
        • Clark K.
        • Vendt B.
        • Smith K.
        • Freymann J.
        • Kirby J.
        • Koppel P.
        • et al.
        The Cancer Imaging Archive (TCIA): maintaining and operating a public information repository.
        J Digit Imaging. 2013; https://doi.org/10.1007/s10278-013-9622-7
        • Hugo G.D.
        • Weiss E.
        • Sleeman W.C.
        • Balik S.
        • Keall P.J.
        • Lu J.
        • et al.
        A longitudinal four-dimensional computed tomography and cone beam computed tomography dataset for image-guided radiation therapy research in lung cancer.
        Med Phys. 2017; 44: 762-771https://doi.org/10.1002/mp.12059
      6. A. Fedorov R. Beichel J. Kalpathy-Cramer J. Finet J.-C.C. Fillion-Robin S. Pujol et al. 3D Slicer as an image computing platform for the Quantitative Imaging Network Magn Reson Imaging 2012;30:1323–41. 10.1016/j.mri.2012.05.001.

      7. Rasband W.S. ImageJ U.S. National Institutes of Health 2018 Bethesda, Maryland, USA. Http://ImagejNihGov/Ij/.

        • Kemerink G.J.
        • Kruize H.H.
        • Lamers R.J.S.
        • van Engelshoven J.M.A.
        Density resolution in quantitative computed tomography of foam and lung.
        Med Phys. 1996; 23: 1697-1708https://doi.org/10.1118/1.597757
        • Ekström M.
        • Sundh J.
        • Schiöler L.
        • Lindberg E.
        • Rosengren A.
        • Bergström G.
        • et al.
        Absolute lung size and the sex difference in breathlessness in the general population.
        PLoS ONE. 2018; 13: 1-12https://doi.org/10.1371/journal.pone.0190876
        • Needham C.D.
        • Rogan M.C.
        • McDonald I.
        Normal Standards for Lung Volumes, Intrapulmonary Gas-mixing, and Maximum Breathing Capacity.
        Thorax. 1954; 9: 313-325https://doi.org/10.1136/thx.9.4.313
        • Obert M.
        • Kampschulte M.
        • Limburg R.
        • Barańczuk S.
        • Krombach G.A.
        Quantitative computed tomography applied to interstitial lung diseases.
        Eur J Radiol. 2018; 100: 99-107https://doi.org/10.1016/j.ejrad.2018.01.018
        • Schroeder J.D.
        • McKenzie A.S.
        • Zach J.A.
        • Wilson C.G.
        • Curran-Everett D.
        • Stinson D.S.
        • et al.
        Relationships Between Airflow Obstruction and Quantitative CT Measurements of Emphysema, Air Trapping, and Airways in Subjects With and Without Chronic Obstructive Pulmonary Disease.
        Am J Roentgenol. 2013; 201: W460-W470https://doi.org/10.2214/AJR.12.10102
      8. Mackin D, Fave X, Zhang L, Fried D, Yang J, Brian Taylor, et al. Measuring computed tomography scanner variability of radiomics features. Invest Radiol 2015;50:757–65. https://doi.org/10.1097/RLI.0000000000000180.

        • Reske A.W.
        • Busse H.
        • Amato M.B.P.
        • Jaekel M.
        • Kahn T.
        • Schwarzkopf P.
        • et al.
        Image reconstruction affects computer tomographic assessment of lung hyperinflation.
        Intensive Care Med. 2008; 34: 2044-2053https://doi.org/10.1007/s00134-008-1175-8
        • Chen-Mayer H.H.
        • Fuld M.K.
        • Hoppel B.
        • Judy P.F.
        • Sieren J.P.
        • Guo J.
        • et al.
        Standardizing CT lung density measure across scanner manufacturers.
        Med Phys. 2017; 44: 974-985https://doi.org/10.1002/mp.12087
      9. Solomon J, Samei E. Quantum noise properties of CT images with anatomical textured backgrounds across reconstruction algorithms: FBP and SAFIRE. Med Phys 2014;41. https://doi.org/10.1118/1.4893497.

        • Nehmeh S.A.
        • Erdi Y.E.
        • Ling C.C.
        • Rosenzweig K.E.
        • Squire O.D.
        • Braban L.E.
        • et al.
        Effect of respiratory gating on reducing lung motion artifacts in PET imaging of lung cancer.
        Med Phys. 2002; https://doi.org/10.1118/1.1448824
        • Liu C.
        • Pierce L.A.
        • Alessio A.M.
        • Kinahan P.E.
        The impact of respiratory motion on tumor quantification and delineation in static PET/CT imaging.
        Phys Med Biol. 2009; https://doi.org/10.1088/0031-9155/54/24/007
        • Rietzel E.
        • Chen G.T.Y.
        • Choi N.C.
        • Willet C.G.
        Four-dimensional image-based treatment planning: Target volume segmentation and dose calculation in the presence of respiratory motion.
        Int J Radiat Oncol Biol Phys. 2005; https://doi.org/10.1016/j.ijrobp.2004.11.037
        • Townsley M.I.
        Structure and composition of pulmonary arteries, capillaries, and veins.
        Compr Physiol. 2012; 2: 675-709https://doi.org/10.1002/cphy.c100081
        • Best A.C.
        • Meng J.
        • Lynch A.M.
        • Bozic C.M.
        • Miller D.
        • Grunwald G.K.
        • et al.
        Idiopathic pulmonary fibrosis: Physiologic tests, quantitative CT indexes, and CT visual scores as predictors of mortality.
        Radiology. 2008; https://doi.org/10.1148/radiol.2463062200