2.1.

Realistic Anthropomorphic Flexible Phantom

Modeling started from the realistic anthropomorphic flexible (RAF) phantom, a full body male phantom developed by Lombardo et al.⁷ using polygonal mesh modeling. This type of geometrical representation is widely used in computer graphic modeling⁸ and describes the phantom with a collection of polygons that share vertices and edges fulfilling certain rules.

The RAF phantom had been validated against the ICRP Publication 110 phantom⁹ and provides a detailed depiction of human anatomy [Fig. 1(a)]. The polygonal mesh format of the phantom allows for further modification of anatomy and posture. This flexibility is important for the accurate modeling of patient position during radiographic examinations, as this can range from an erect chest posterior anterior (PA) to a bedside anterior posterior examination. For current chest imaging applications, only the organs in the thorax region of the RAF phantom have been developed further. Additional modifications included the modeling of a more detailed lung background¹⁰ with bronchial trees, pulmonary arteries, and pulmonary veins [Figs. 1(b) and 1(c)].

Fig. 1

(a) Picture of the RAF mesh phantom; (b) mesh version of the RAF phantom adapted for chest imaging applications; (c) voxelized version of the RAF phantom; and (d) different versions of the RAF phantom representing (from left to right) the standard male, the modified overweight male, and a female.

The RAF phantom represents a human male with a body mass index (BMI) of $24\mathrm{kg}/{\mathrm{m}}^2$. Additional versions of the phantom were created to represent an overweight ($\mathrm{BMI}=29\mathrm{kg}/{\mathrm{m}}^2$) and obese ($\mathrm{BMI}=40\mathrm{kg}/{\mathrm{m}}^2$) male as well as a female version with a BMI of $29\mathrm{kg}/{\mathrm{m}}^2$ [Fig. 1(d)]. The different body types were created by modifying the external shape of the phantom. Graphic modeling software 3ds Max (Autodesk, USA) was used for this purpose. The soft selection tool was employed as this allows the vertices of the meshes to be deformed in an organic way. Figure 1(d) shows an example of the soft selection; the red color represents the meshes explicitly selected, whereas the remaining colors represent the vicinity meshes. As the selected meshes are transformed (i.e., translated, rotated, and/or scaled), the elements in the vicinity are drawn along smoothly, and this effect decreases with distance or the “strength” of the selection. The internal organs of the phantom were kept unchanged. A study by Lemanowicz et al.¹¹ found the level of patient obesity to have the highest correlation with the chest soft tissue thickness; thus, it was considered correct to keep the internal organs unchanged for the first set of models.

2.2.

Modeling of the Pathologies from CT Scans

Ethical approval was requested for a retrospective study making use of patient CT images; the ethical committee waved any dedicated patient consent. A 3D modeling methodology was implemented to create the COVID-19 disease within the RAF phantom. CT images of patients suspected of COVID-19 confirmed by CT scans were used as reference for the development of the pathology models. The CT scans had been acquired with a low-dose thorax protocol using 120 kVp, 1.2 mm pitch, and 46 mAs (tube current modulation off). CT image voxel size ranged from 0.68 to 1.00 mm and slice thickness was 3 mm.

A range of cases was selected by a radiologist to ensure that the different stages of the disease would be covered. Areas of ground glass opacities (GGO) and consolidation were segmented manually using ImageJ¹² [Fig. 2(a)]. The segmentation was carried out by a medical physicist but was closely guided by a radiologist. The segmented pathology was then converted to a binary stack image in which the pathology was colored white and the background black. Next, the marching cubes algorithm¹³ was applied to extract polygonal meshes from the isosurfaces of the three-dimensional pathological structures [Fig. 2(b)]. This conversion from CT voxels to meshes was done using the 3D visualization library in ImageJ.¹⁴ A resampling factor of 1 was chosen for the marching cubes, so the meshes would replicate the voxel structures without a loss of resolution. The mesh volume of the pathology was then exported as an OBJ file containing the coordinates of the vertices in 3D space; this in turn defined the shape and size of the surface of the segmented pathology.

Fig. 2

Workflow of the methodology followed to model the pathologies: (a) patient CT slice with segmented pathology (green); (b) 3D surface of the pathology converted to mesh format; (c) slice of the 3D mesh model of the pathology (green) inside RAF phantom lungs; and (d) slice of voxelized RAF phantom with highlighted pathology (green).

The OBJ lesion model was imported into 3ds Max, where two steps were used to correct the meshes and generate consistent models. First, the voxel scaling effect due to the segmentation process in anisotropic CT data was smoothed. Artifacts generated by the marching cubes algorithm were removed using the TurboSmooth and Relax modifiers. The latter were sparingly applied to avoid loss in the model volume. Second, the shape of the pathology was modified to fit the RAF phantom lungs, while respecting the volume of the segmented disease [Fig. 2(c)]. This step was carried out using the free form deformation modifier. The deformations applied in this step required manual intervention since the shape of the lung and its pathology changes from case to case. Non-isotropic scaling was utilized. The original spatial distribution of the pathology over the lung was preserved, and changes in the volume ratio disease/lung ($V_{\text{segmented}\_\text{pathology}}/V_{\text{patient}\_\text{lung}}$) were kept below 15%. The set of models represented the typical distributions of the disease: predominantly in the lower lobe, multifocal, peripheral, and bilateral.¹⁵ To quantitatively assess differences between the segmented lesion and the final model used in the phantom, the Hausdorff distance (HD) was calculated for the mesh models. Meshlab software,¹⁶ which samples a set of points over the mesh and finds the closest point in the reference mesh, was used.

The modeling methodology enables the representation of different grades of severity of lung involvement by changing the x-ray attenuation in the simulated pathology or by changing its size and distribution. As a proof of concept for the development of different lesion severities from the reference set, the mesh of one of the segmented lesions was modified by changing its size and shape.

An example of a model in which the x-ray attenuation of the pathology was changed to study the impact on lesion detection was also simulated. Finally, the influence of BMI on lesion detection is also illustrated for one disease model simulated in the different BMI and gender realizations of the phantom. These examples are shown in the results (Sec. 3.5).

2.3.

Creation of Voxelized Phantoms

The polygonal mesh models were then exported separately and loaded into voxelization software.⁷ The latter software tool is based on the work of Laine¹⁷ and is optimized such that the mass and thickness of the organs is preserved. The algorithm uses a conservative eight-separating voxelization method by calculating intersections between the triangles and quads from the meshes with a cube used as the test object.¹⁷ The voxel resolution used for the models was $0.5\times 0.5\times 0.5{\mathrm{mm}}^3$, which was a compromise that maintained fine detail within the phantom anatomy and yet kept the computational load reasonably low. Once all of the organs and pathologies were voxelized; they were combined and assigned an ID number using a script developed for ImageJ.¹² The process started from the thoracic wall and went in an increasing hierarchy to the bones, diaphragm, lungs, heart, up to the pulmonary veins, vessels, and the bronchi (highest priority) [see Fig. 2(d)]: IDs were overwritten where there was superposition. The GGO ID was added after the lungs since GGO represent an area of increased opacity in which bronchial structures and vessels are still visible. The consolidation ID was added last, as consolidations obscure all pulmonary structures within some region¹⁸ (see Fig. 3).

Fig. 3

Slice of mesh model of the lung containing GGO and consolidation regions. GGO are regions of increased opacity in the lungs where pulmonary structures are still visible, whereas in the consolidation the pulmonary structures are obscured.

Finally, the volume ratio of the segmented pathology to that of the patient lung ($V_{\text{segmented}\_\text{pathology}}/V_{\text{patient}\_\text{lung}}$) was compared with the volume ratio of the voxelized pathology model to that of the phantom lung ($V_{\text{final}\_\text{pathology}}/V_{\text{phantom}\_\text{lung}}$). This served as a validation of the method to ensure that the volume ratio of pathology/lung was kept similar, regardless of differences between the patients and the RAF chest phantom.

2.4.

Generating Radiographic Images Using a Simulation Framework

An imaging chain simulation platform developed in-house was used to generate radiographic images of the models. The simulation platform includes all elements from the imaging chain: x-ray source, patient model, antiscatter grid,¹⁹ and detector. The simulations combined a ray tracing algorithm²⁰ and Monte Carlo simulations using the PENELOPE/penEasy^21,22 transport code. The ray tracing was used to create noise free, high-resolution primary projections (i.e., images generated by unscattered photons) in relatively short computation times. PenEasy was used to create scatter images by tracking photons that underwent energy loss or angular deflection before arriving in the detector. The random nature of the Monte Carlo simulations required long computational times to obtain low uncertainty values. To reduce the computational time, pixels of $5\times 5{\mathrm{mm}}^2$ were used in the scatter image. This was considered sufficient to represent the variations of the scatter radiation since these images do not contain fine details.²³ The uncertainty of the Monte Carlo calculations was kept below 3%. The materials and proportions of the mixtures used for the different phantom organs were taken from ICRP Publications 89²⁴ and 110.⁹ GGO and consolidations are areas of increased levels of x-ray attenuation due to the presence of fluid in the lungs. To obtain realistic pathology densities, the Hounsfield unit (HU) histograms of the corresponding segmented lesions in the CT images were measured and subsequently converted to density ($\mathrm{g}/{\mathrm{cm}}^3$) using the data of Schneider.²⁵ A random distribution of densities was used in the pathology region of the phantom. A random value generator was used to assign different attenuation coefficients to each voxel in the pathology region: each density value was given the probability of the corresponding HU in the CT images. To achieve this, the probability distribution was taken from the normalized HU histogram of the segmented lesion.

Realistic levels of sharpness and noise, relevant to a typical x-ray detector for chest projection imaging, were then applied to the simulated phantom images. Sharpness and noise were quantified using the presampling modulation transfer function (MTF) and the normalized noise power spectrum (NNPS), respectively. These metrics were measured from a Carestream flat panel CsI digital detector with a beam quality of 120 kVp and 9 cm PMMA (i.e., lung equivalent thickness¹⁰) at the tube exit. The presampling MTF was measured using a version of the edge technique described by Samei et al.²⁶ A square tantalum edge test object of dimension $50\times 50{\mathrm{mm}}^2$ and 1-mm-thickness was placed at the center of the image receptor and oriented manually to obtain an angle of $\sim 3\mathrm{deg}$ between the edge and the detector matrix. The MTF images were acquired at a detector air kerma (DAK) of $8\mu \mathrm{Gy}$. The NNPS was measured from flat field images acquired at different dose levels, with DAK values ranging from 0.7 to $25\mu \mathrm{Gy}$.²⁷ A region of $1024\times 1024\text{pixels}$ was extracted from the image center and a 2D polynomial fitted to and then subtracted from this region to reduce the influence of large area non-uniformities on the final NNPS.²⁸ Half-overlapping regions of interests of dimension $128\times 128\text{pixels}$ were then extracted from this region, from which the 2D NNPS was calculated.²⁹

The sharpness and noise, as characterized by the MTF and NNPS, respectively, were then applied as follows. The 1D presampling MTF was fitted and expanded to a 2D MTF using weighting matrices described in the literature.³⁰ The fast Fourier transform (FFT) of the simulated phantom image was then multiplied with the 2D MTF, and an inverse FFT applied to obtain a noise-free image with a realistic level of blurring.²³ Image noise was then applied to the blurred images via the NPS following the methods of Båth et al.³¹ and Mackenzie et al.³² The method consists of forming a noise image using the different components of the NNPS, i.e., electronic, quantum, and fixed pattern noise. Three random Gaussian white noise images with zero mean and unit standard deviation were weighted by these components. These images were individually multiplied with the blurred phantom images and added to form the total noise image. The final simulated x-ray projection was obtained by summing the noise image and the blurred image. Radiographic images of all of the models were generated at 120 kVp, grid in, and 180 cm source to detector distance, settings commonly used in thorax PA examinations. A DAK level of $\sim 8\mu \mathrm{Gy}$, which is higher than in the clinical protocol, was simulated. A PA exposure setting with grid and a higher dose level were selected to obtain high-quality images, in which the pathology was clearly visualized. Finally, clinical image processing corresponding to an adult thorax examination was applied to the images using image processing software MUSICA (Agfa, Belgium).

The imaging chain simulation was validated by comparing the signal difference to noise ratio (SDNR) values from real and simulated images of a test object consisting of a PMMA block (10 and 20 cm thickness) and an aluminum detail (2 mm thickness). Relative differences between SDNR in the real and simulated images were below 10%. This validation was performed over a range of tube voltages and dose levels.

2.5.

Assessment of Task Realism

To assess the realism of the pathologies in the RAF phantom, the simulated radiographic images were presented to three thorax radiologists, trained during the outbreak to diagnose COVID-19 suspected patients. A reader study in which the observers were asked to score nine different images for each of the COVID-19 models inserted within the standard BMI male RAF was set up. The score was carried out according to three realism criteria: question 1: realism of the lung background, question 2: realism of the lesions in terms of appearance, and question 3: realism of the lesions in terms of position within the lungs. A five-point scale was used for all criteria: (1) not at all realistic: critical elements that affect the realism, (2) poor: obvious elements that may affect the realism, (3) adequate: minor elements that did not affect the realism, (4) good: minimal unrealistic elements, and (5) very realistic: no unrealistic elements. In addition, the readers were asked to write a description of the pathology that was seen. Images were displayed using Viewdex software.³³ Contrast and brightness levels could be adjusted if desired, and no time limit was imposed.

As an additional validation of the models, the images were uploaded to the AI software Lunit INSIGHT CXR for COVID-19 (Lunit, South Korea).¹⁸ This software can detect areas of consolidation or GGO in the chest and is intended to support the interpretation of CXR of suspected COVID-19 cases. The software analyses the images and reports an abnormality score given by the likelihood of the presence of the detected lesion, i.e., low (0% to 15%), moderate (16% to 50%), and high (51% to 100%).

3.1.

Pathology Models

Nine COVID-19 disease models were created, covering the typical manifestations of lung involvement distribution characteristics of COVID-19 pneumonia. Eight of the models were created from the segmented lesions in CT images and one extra case was created by modifying the mesh of one of the segmented models. The HU histograms obtained for the lesions of the eight segmented CT datasets are shown in Fig. 4. The differences in HU are related to the stage of the disease. Although GGOs lie in the range $-800$ to $-400$, higher opacities like consolidations reach 100 HU, consistent with data reported by Lanza et al.³⁴

Fig. 4

HU distribution within the lesions for the pathologies modeled from the segmentation of CT images (cases 1 to 8).

Figure 5 shows a comparison of different slices of the CT data [(a)–(c)] to the corresponding slice of the voxelized phantom [(d)–(f)] for case 4. The pathology voxels are highlighted in green in all of the images.

Fig. 5

(a)–(c) Comparison of CT slices of a patient (case 4) and (d)–(f) the respective slices of the voxelized RAF phantom. Voxels corresponding to the pathology are highlighted in green.

Different stages of the disease were modeled, from the more subtle in the initial stage to more advanced. Figures 6(a)–6(c) show the mesh models of the RAF phantom lungs including the pathologies from cases 6, 4, and 7, respectively. As illustrated, the level of lung involvement changes from case to case; the ratios of pathology volume to lung volume of the voxelized versions of these models can be found in Table 1.

Fig. 6

3D mesh models of the pathologies (green) modeled within the lungs of the RAF phantom. Different levels of lung involvement can be seen, from subtle (a) case 6 to more prominent in (b) case 4 and (c) case 7.

Table 1

Volume ratios of pathology to lung for the segmented CT images and those for the corresponding phantom models. Relative deviation between the developed models and the patient data is also shown. Mean Hausdorff distance calculation for the mesh models of cases 1 to 8. No CT data are available for case 9 since it is a modified version of case 4.

3.2.

Mesh and Volume Comparison

Table 1 shows the volume ratios of pathology/lungs for the developed models. As can be seen, the percentage of lung involvement ranges from 2.2% to 38.2%. The pathology/volume ratios for the segmented CT images and the relative differences compared with the voxelized models are also shown. As observed, maximum differences between the initially segmented lesion and the corresponding final phantom version stayed below 16%. Volume ratios were usually smaller when the pathology was placed in the phantom. No comparison is shown for case 9 since this model is a modification of case 4; thus, no reference volume ratio is available.

The mean HD between the original mesh and the mesh adapted to the RAF phantom can be found in Table 1. Mean HD values were below 0.64 cm.

3.3.

Simulated Radiographic Images

Figures 7(a)–7(i) show the set of simulated radiographic images of the RAF phantom featuring the nine pathology models. For each case, a projection is also shown with the voxels containing the pathology highlighted in green. As observed, the pathologies in all cases are bilateral, often located in the periphery of the lungs and involving several lobes. The degree of spread of the disease is easily noticeable in each of the cases.

Fig. 7

Simulated images of the RAF phantom including the COVID-19 models created. From (a)–(i) cases 1 to 9, respectively. Below each image, an equivalent version with pathology regions highlighted in green is shown.

3.4.

Assessment of Task Realism

Table 2 presents the percentage of cases marked as at least adequate or at least good for the individual observer and the mean realism criteria of all three observers used in the reader study. For questions 1, 2, and 3 (Q1, Q2, and Q3, respectively), the average percentages of cases marked at least adequate were 100%, 92%, and 96%, respectively. The average percentage of cases marked at least good were 59%, 54%, and 65% for Q1, Q2, and Q3, respectively. The mean value for the three quality criteria was at least adequate for 96% of the cases and at least good for 50% of the cases.

Table 2

Percentage of cases rated by radiologists as at least adequate and at least good for each of the realism criteria in the images. Q1: realism of lung background, Q2: realism of lesion (appearance), and Q3: realism of lesion (position).

For readers 1, 2, and 3, 92%, 96%, and 100% of all scores given to the images, respectively, were at least adequate; they detected only minor unrealistic elements that did not affect the general realism of the models. Moreover, 67%, 52%, and 59% of all scores were at least good for readers 1, 2, and 3, respectively. A moderate agreement was found between the readers with an intraclass correlation coefficient (ICC) of 0.5 (95% $\mathrm{CI}=0.10$ to 0.84).

Radiologists 1, 2, and 3 detected 78%, 89%, and 89% of all COVID-19 pathologies, respectively. Case 6 was missed by all three readers, and case 5 was missed by reader 1. Although reader 3 had marked an area of opacities in this lung, he was uncertain about its presence. These two cases, in fact, represent the more subtle simulated pathologies. Further investigation of case 5 revealed, that after three PCR tests, this patient was negative for COVID-19, and it represents a rare case where CT and PCR arrived at different conclusions. This highlights the difficulty in the radiological practice where lesions with spatial distributions typical for COVID-19 and HU distributions like COVID-19 cases can in fact be negative.

The AI algorithm was able to identify consolidation areas in 56% of the cases. The algorithm successfully identified pathologic areas in cases 3, 4, 7, 8, and 9. For the rest of the cases, no lesions were detected. Figures 8(a)–8(e) show simulated images with the pathology highlighted in green and the corresponding output from the AI software for cases 3, 4, 7, 8, and 9, respectively. The heat map displayed from the AI software output represents the likelihood of the presence of the lesion as detected by the software. The abnormality scores reported by the AI were 90%, 93%, 77%, 92%, and 92% for cases 3, 4, 7, 8, and 9, respectively.

Fig. 8

Simulated images for cases (a) 3, (b) 4, (c) 7, (d) 8, and (e) 9. In the upper part is the image with the pathology highlighted in green and below is the corresponding result from the AI software. The heat maps in the bottom images represent the likelihood of the detected lesion to be suggestive for COVID-19.

3.5.

BMI and Pathology Modifications

Figure 9 presents the images from the RAF phantom with different body types, namely (a) female, (b) overweight male, and (c) obese male. The same pathology, taken from case 4, was used for all of these images [Fig. 7(d)]. As observed, the visualization of the pathology is affected by the increased size of the phantom and by the presence of breasts. These images were analyzed by the AI COVID-19 detection software, and the corresponding output is shown below each case. The different body types had a clear influence on the detection of the pathology, as can be seen in the heat maps of the AI software. The abnormality scores were 88%, 78%, and 56% for the female, male overweight, and male obese, respectively, compared with 93% obtained for the same case in the standard BMI male.

Fig. 9

Images of the RAF phantom (with the pathology of case 4) with different body types (a) female ($\mathrm{BMI}=29$); (b) male overweight ($\mathrm{BMI}=29$); and (c) male obese ($\mathrm{BMI}=40$); below each image, the respective output from the AI software is displayed.

Figure 10 shows the x-ray images of two versions of case 8: (a) represents the original pathology and (b) is a modified version where lower attenuation in the lungs was simulated to represent a more subtle case. These images were uploaded to the AI software for analysis. The respective results are displayed below each simulated image. Heat map intensity falls for the more subtle case, with the abnormality score going from 92% in the original model to 77% in the modified subtler version.

Fig. 10

Example of different stages of the disease for case 8. (a) The original models and (b) a more subtle stage with less attenuation in the lungs. The corresponding AI output is shown.