WO2022164374A1

WO2022164374A1 - Automated measurement of morphometric and geometric parameters of large vessels in computed tomography pulmonary angiography

Info

Publication number: WO2022164374A1
Application number: PCT/SE2022/050089
Authority: WO
Inventors: Ali Teymur KAHRAMAN; Tobias SJÖBLOM; Tomas FRÖDING
Original assignee: Kahraman Ali Teymur; Sjoeblom Tobias; Froeding Tomas
Priority date: 2021-02-01
Filing date: 2022-01-28
Publication date: 2022-08-04

Abstract

The present disclosure relates to a computer-aided detection, CADe, method performed in a computing device for locating and measuring large vessels in a subject undergoing a computed tomography, CT, scan, comprising obtaining CT data from a consecutive CT scan of thorax of the subject, the CT data comprising an image stack of slices forming a CT volume of images, preparing the CT data for segmentation tasks, locating carina trachea in the CT volume, detecting descending aorta at carina level, detecting ascending aorta, detecting pulmonary trunk, and detecting an apical level of pulmonary valve. The method further relate to a computing device for performing the method, and a computer-aided detection system comprising the computing device.

Description

AUTOMATED MEASUREMENT OF MORPHOMETRIC AND GEOMETRIC PARAMETERS OF LARGE VESSELS IN COMPUTED TOMOGRAPHY PULMONARY ANGIOGRAPHY

TECHNICAL FIELD

The present disclosure relates automatic measurements of morphometric and geometric parameters of large vessels. More specifically, the proposed technique relates to methods for automated measurement of large vessels in computed tomography pulmonary angiography, for detecting or diagnosing conditions or disorders present in the subject undergoing the scan. The disclosure further relates to a computing device for performing the methods, and a computer-aided detection system comprising the computing device. The disclosure also relates to computer programs and carriers thereof.

BACKGROUND

Over the last few decades, the technological development of computer systems and medical imaging systems has increased the quality of medical images and has decreased the image acquisition time notably. This technological advance has made the use of medical images indispensable in the detection and diagnosis of various diseases. On the other hand, this rapid development has increased the number of cases (medical images) to be examined by radiologists, and the workload of radiologists has therefore increased. Accordingly, the multitude of medical images obtained makes it challenging for radiologists to examine and report cases on time.

Performing medical imaging, such as computed tomography (CT) imaging, of the chest of patients is a well-established method for examining patients in need thereof. These examinations may also contain more information than the examining radiologist can interpret primarily because of time constraints. The radiologist seeks signs of specific conditions, and may therefore overlook other information present in the images. For example, enlargement of the aorta or pulmonary trunk may be indicative of disorders such as aneurysm or pulmonary embolism. Detecting such abnormalities may also be challenging for the examining radiologist, or they may be down-prioritized due to the radiologist seeking signs of other disorders. Accordingly, diagnosis of severe conditions may be delayed. SUMMARY

An object of the present disclosure is to provide methods and devices which seek to mitigate, alleviate, or eliminate the above-identified deficiencies in the art and disadvantages singly or in any combination. This object is obtained by computer-aided detection, CADe, method performed in a computing device for locating and measuring large vessels in a subject undergoing a computed tomography, CT, scan , the method comprising: obtaining CT data from a consecutive CT scan of thorax of the subject, the CT data comprising an image stack of slices forming a CT volume of images, preparing the CT data for segmentation tasks, locating carina trachea in the CT volume, detecting, descending aorta at carina level, detecting ascending aorta, detecting pulmonary trunk, and detecting an apical level of pulmonary valve.

In some aspects the method further comprises determining a representative diameter of the ascending aorta by segmenting the ascending aorta, tracking the ascending aorta slice by slice from the first located slice to the carina trachea, wherein a diameter of the segmented ascending aorta is measured in every slice, and calculating a representative diameter of the ascending aorta as a mean of the measured diameters from each slice.

In some aspects the method further comprises determining a representative diameter of the pulmonary trunk by segmenting the pulmonary trunk, and tracking the segmented pulmonary trunk, slice by slice, from carina level to pulmonary valve and back to carina level, wherein a diameter of the pulmonary trunk is measured in each slice, and calculating a representative diameter of the pulmonary trunk as a mean of the measured diameters from each slice.

In a further aspect is provided a computer-aided detection, CADe, method performed in a computing device for locating and measuring large vessels and detecting ascending aorta hypertension and/or pulmonary trunk hypertension in a subject, the method comprising: performing any of the method above or below to determine a representative diameter of an ascending aorta and/or pulmonary trunk in a subject, comparing the determined representative diameter(s) with a preset threshold, determining that ascending aorta hypertension and/or pulmonary trunk hypertension is present in a subject if the determined respective representative diameter(s) is greater than the preset threshold.

In a further aspect is provided a computing device comprising a memory for storing instructions, and processing circuitry for executing the instructions, wherein the processing circuitry is configured to perform the methods of the invention.

In a further aspect is provided a computer aided detection system comprising a computing device, configured for performing a computer-aided detection, CADe, method of locating and measuring large vessels in a subject undergoing a computed tomography, CT, scan, the computing device comprising, a memory for storing instructions and processing circuitry configured to cause the computing device: to obtain CT data from a consecutive CT scan of thorax of the subject, the CT data comprising an image stack forming a CT volume of images, to prepare the CT data for segmentation tasks, to locate carina trachea in the CT volume, to detect descending aorta at carina level, to detect ascending aorta, to detect pulmonary trunk, and to detect an apical level of pulmonary valve. The computer-aided detection system of is further configured to perform the methods of the invention.

According to some aspects, the computer-aided detection system is configured for determining a representative diameter of a pulmonary trunk in a subject, the system comprising: a processor configured to determine a representative diameter of a pulmonary trunk in a CT volume of images from a CT scan of the subject by tracking a segmented pulmonary trunk from carina level to pulmonary valve and back to carina level, wherein a diameter of the pulmonary trunk is measured in each slice, and a representative diameter of the pulmonary trunk is calculated as a mean of the measured diameters from each slice.

According to some aspects, is provided a computer program comprising computer program code which, when executed in a computing device, causes the computing device to execute the methods according to the invention, and a carrier containing the computer program, wherein the carrier is one of an electronic signal, optical signal, radio signal, or computer readable storage medium. BRIEF DESCRIPTION OF THE DRAWINGS

Figure 1 shows an illustration of the CT scanner principle.

Figure 2 illustrates a voxel in the CT slice.

Figure 3 is an illustration an image acquisition method of the present invention.

Figure 4 is a flowchart of an exemplary process for a CADe method performed in a computing device for locating and measuring large vessels in a subject undergoing a CT scan.

Figure 5 illustrates finding the orientation of the subject by using the angle between the major axis of the subject's image and the x-axis of the image plane.

Figure 6 shows the algorithmic steps of a first patient orientation calculation method.

Figure 7 A and B show the algorithmic steps of a second patient orientation calculation method.

Figure 8 illustrates finding an accurate orientation of the CT scan.

Figure 9 shows artificial rays defining search space 1 (the left image) and search space 2 (the right image).

Figure 10 illustrates the location of the aortic arch.

Figure 11 illustrates the search space for the pulmonary trunk (PT) in an image slice.

Figure 12 illustrates a tracking area for the pulmonary trunk.

Figure 13 shows applying the Hough transform to a segmented component of a pulmonary artery for locating the apical level of pulmonary valve.

Figure 14 shows a block diagram of an example computing device of the invention.

Figure 15 is a demonstration of manual measurements by radiologist.

Figure 16 illustrates examples of CT pulmonary angiography image qualities.

Figure 17 illustrates a flowchart of the CADe algorithm using images.

Figure 18 shows an example of successful segmentation of the trachea.

Figure 19 is an example of successful segmentation of the descending aorta.

Figure 20 is an example of successful segmentation of the ascending aorta.

Figure 21 is an example of successful segmentation of the pulmonary trunk.

DETAILED DESCRIPTION Aspects of the present disclosure will be described more fully hereinafter with reference to the accompanying drawings. The computing device and methods disclosed herein can, however, be realized in many different forms and should not be construed as being limited to the aspects set forth herein. Like numbers in the drawings refer to like elements throughout.

The terminology used herein is for the purpose of describing particular aspects of the disclosure only, and is not intended to limit the disclosure. As used herein, the singular forms "a", "an" and "the" are intended to include the plural forms as well, unless the context clearly indicates otherwise.

In some embodiments the non-limiting term "medical imaging" is used. Medical imaging techniques comprise imaging techniques for clinical applications, such as X-ray, Computed Tomography (CT), Magnetic Resonance Imaging (MRI), Ultrasound, Nuclear Medicine (Single Photon Emission Computed Tomography (SPECT), Positron Emission Tomography (PET)) and Endoscopy.

The term "computed tomography" or "CT" as used herein may refer to any type of X-ray computed tomography (X-ray CT) or computerized axial tomography scan (CAT scan), and may be performed by several different systems, such as Standard (Conventional) CT, Helical (Spiral) CT (HCT), High-Resolution CT (HRCT), Multiple-row Detector CT (MDCT), and Dual Source CT (DSCT).

A "CT scan" of a subject as used herein may refer to performing all steps necessary to achieve a volume of images from the subject, including placing the subject on a motorized table, moving the table into gantry while X-ray sources rotate around the table. Recoding the subject's body snapshots, which are exited from X-ray beams, using detectors, and sending the snapshots to a computer to reconstruct them into axial slices (image slices), where each snapshot is turned into an image matrix (slice), which is made of numerical values. A "CT scan" may either refer to a single scan (snapshot) performed by the CT on a single slice of the body of the subject, or to a complete CT scan of the subject's body part of interest comprising a plurality of consecutive images. The CT scan give rise to an image stack forming a CT volume of images. The CT volume of images may also be referred to as "CT data". The "CT data" may also be referred to as comprising a CT volume of images or CT image volume.

The terms "image", "slice" and "image slice", as in for example "CT image" or "CT slice", are used interchangeably herein and refers to a single image of a single slice of the body captured during the CT scan, a single two-dimensional (2D) image plane. Accordingly, an "image stack" consists of a plurality of consecutive CT scans comprising a plurality of images, also referred to as CT "slices" or CT "volume images", since the image stack of plurality of consecutive images together span a volume of the scanned subject.

Segmenting and extracting may be used interchangeably herein and refers to grouping a set of pixels in one segment. A 3D segmentation of an object done in one direction. Once a region of interest is segmented (or extracted) in a CT slice, checking if this segmented region exist in next or previous CT slice (depending on which direction one needs to go, cranial to caudal or caudal to cranial) is a tracking procedure. At the end of the tracking procedure, are the volumes of interest (areas adjacent to each other and sequential on the z-axis).

In some instances, a distance or area is referred to as a number of pixels. "Pixel Resolutions" (width x height) is the number of total pixels in the image represented as width x height (or rows x columns, or matrix size). The pixel resolutions in the cohort of the present disclosure are 512 x 512 (meaning 512 pixels horizontally by 512 vertically). The standard resolution used today by the CT manufacturer is 512 x 512 and can be found at the rows and columns meta tag of the CT examination. "Pixel Spacing" is the distance between two neighbor pixels in mm.

The pixel spacing value can be found at the meta tag (pixel spacing tag) of the CTPA examination. There is no standard pixel spacing value. Therefore, the pixel spacing values can be different regarding the CT device manufacturer. In our cohort (n=700), pixel spacing varies between 0.58 mm to 0.9 mm (mean =0.72 mm, median = 0.70 mm, mode = 0.70 mm). Thus, a distance or an area defined by a number of pixels may be translated into metric scale, the formula of the pixel to the metric scale being = (total number of pixels) x square of pixel spacing value. Thus, to convert pixels numbers to metric scale, one only needs pixel spacing value (a.k.a. pixel size). For example, if the values in the pixel spacing tag are 0.70/0.70, it corresponds to a distance between the two neighbor pixels of 0.70 mm. Thus, if the total number of pixels is 5 and the pixel spacing value is 0.7 mm then, 5 pixels area equals = 5 x (0.7 x 0.7) = 2.45 mm².

Where a distance or area is referred to as a number of pixels, it is understood in this regard that the distance or area in number of pixels apply when the same resolution is used as in the examples of the present disclosure. Thus, if a different resolution is used, the distance or area must be translated. In some instances, a distance of 11 or 15 pixels are defined. In other instances, 5, 10, 15, 100, 350, 1200, and 1750 pixels are mentioned. The number of pixels in these regards may be translated into the metric system in two ways.

In a first way, since pixel spacing values vary, we need to use a constant spacing mean value (0.72 mm). Square of pixel spacing value = 0.72 x 0.72 mm ~ 0.5 mm² . Thus, 5 pixel -> 5 x 0.5 =2.5 mm², 10 pixel -> 10 x 0.5 = 5 mm², 10 pixel -> 11 x 0.5 = 5,5 mm², 15 pixel -> 15 x 0.5 = 7.5 mm², 100 pixel -> 100 x 0.5 = 50 mm², 350 pixel -> 350 x 0.5 = 175 mm², 1200 pixel -> 1200 x 0.5 = 600 mm², and 1750 pixel -> 1750 x 0.5 = 875 mm². The first way is the easiest one, but problematic for the extreme pixel spacing values (ex.: 0.58 mm and above 0.8 mm). Thus, as a second option, a conversion table may be provided in translating pixel numbers to metric scale, see Table 1.

Table 1

Medical imaging is a well-established method of examining subjects in need thereof, i.e. patients suspected to suffer from some kind of condition or disorder. Rapid development of computer systems and medical imaging systems has increased the number of medical images to be examined by the clinicians, and the workload of clinicians has grown rapidly. The multitude of medical images obtained makes it difficult to examine and report cases by clinicians. Computer-aided systems (also called computer-assisted systems) can be used to overcome this problem. Computer-aided systems help clinicians and reduce their workload and increase their efficiency. Advanced computer-aided systems in medical imaging can be divided into three phases, an interpretation phase, a diagnosis phase, a prognosis phase. The interpretation phase contains several image processing and image analysis tasks to segment medical images into the meaningful components, e.g. answering the question of "are these lung fields segmented from CT volume image correctly and are there any abnormalities existing in the segmented component?" The interpretation phase of computer-aided systems is also known as computer-aided detection (CADe) system. In the diagnosis phase, the diseases on the pre-segmented component are diagnosed, e.g. answering the question of "is this tissue, which is extracted from pre-segmented lung fields, a lung cancer?" This phase is also known as computer-aided diagnosis (CADx) system. In the prognosis phase likely outcome of the diagnosed disease is calculated, e.g. answering the question of "what is the survival rate of patients with this kind of lung cancer?" This phase is known as computer-aided prognosis (CAP) system, and CAP is a newer system than the CADe and CADx systems.

The present disclosure will address both the interpretation phase and the diagnosis phase, by automatically segmenting and extracting major parts of the thoracic cavity from chest CT volume images of a subject, and evaluating these to detect abnormalities which may be used as decision support for diagnosis of several conditions.

Pulmonary embolism (PE) is a blockage of an artery in the lungs by a substance that has moved from elsewhere in the body through the bloodstream (embolism). PE is a severe disease which threatens the public health and is associated with high mortality and morbidity rates. Computed tomography pulmonary angiogram (CTPA) is a medical diagnostic test that employs an imaging method, computed tomography angiography, to obtain an image of the pulmonary arteries, and which is an ultimate gold standard for clinical diagnosis of PE. However, such a scan is often only performed if PE is suspected clinically.

Computed Tomography Pulmonary Angiography (CTPA) is commonly used for diagnosis of PE. CTPA operates on computer tomography with iodine-containing contrast agent. The contrast agent increases the density of pulmonary arteries and makes them appear bright. On the contrary, it makes blood clots appear dark in the pulmonary arteries. Thus, PE can be recognized by filling defects in CTPA images, lodine-containing contrast agent increases the density of blood to enhance the appearance of vascular structures in CT scans. X-rays are taken by CT device after the iodine-containing contrast agent is injected into the patient artery. The patient receives an intravenous injection of an iodine-containing contrast agent at a high-rate using an injector pump. Images are acquired with the maximum intensity of radio-opaque contrast in the pulmonary arteries. This can be done using bolus tracking, which is a technique to optimize timing of the imaging. A small bolus of radio-opaque contrast media is injected into a patient via a peripheral intravenous cannula. Depending on the vessel being imaged, the volume of contrast is tracked using a region of interest (ROI) at a certain level and then followed by the CT scanner once it reaches this level. Images are acquired at a rate as fast as the contrast moving through the blood vessels. This method of imaging is used primarily to produce images of arteries, such as the aorta, pulmonary artery, cerebral, carotid and hepatic arteries.

Pulmonary hypertension (PH) is a type of high blood pressure that affects the arteries in your lungs and the right side of your heart. Identification of pulmonary artery dilation on CT scans of the chest of a patient may predict pulmonary hypertension. Pulmonary hypertension may give rise to an increased size of the pulmonary trunk, which may be indicative of a PE in the subject. Accordingly, measuring the average diameter of the pulmonary trunk in a subject undergoing a medical imaging procedure could be used as a diagnostic tool for PE. However, measuring the diameter of the pulmonary trunk manually by the clinician is difficult and labor intensive, and is only made when indications of abnormalities are known. To scan all images of subjects undergoing medical imaging of the chest/thorax region for various reasons also for other abnormalities, such as increased diameters of the pulmonary trunk or aorta, would thus not be possible manually. Accordingly, many theoretically possible diagnostic evaluations of the obtained images are not performed. Accordingly, measurements of morphometric and geometric parameters of mediastinal vessel structures in CTPA examinations can aid diagnosis of cardiovascular conditions like pulmonary embolism and pulmonary hypertension. However, manually measuring these features is time-consuming and a fully automated computer aided system is warranted. Computed tomography CTPA is a medical imaging procedure where an intravenous contrast agent is used to enhance visualization of the thoracic blood vessels, in particular the pulmonary arteries. CTPA is the current reference imaging method for diagnosis of PE. In the PE CTPA reporting process, depending on the clinical situation and local tradition, radiologists measure diameters, volumes, and contrast densities of anatomical structures for the purpose of quality control, diagnosis and risk stratification. Measures such as standard deviation (SD) of CT numbers in a given region of interest (ROI) can be used for image quality assessment. Some of the morphometric and geometric parameters can be of clinical significance aiding diagnosis of conditions such as PE and PH. However, manual measurement of mediastinal structures is time consuming.

Several solutions have been proposed to segment and measure mediastinal structures using CTPA examinations to overcome this problem (see e.g. Saur SC, Kuhhnel C, Boskamp T, Szekely G, Cattin P. Automatic ascending aorta detection in CTA datasets. In: Bildverarbeitung fur die Medizin 2008. Springer; 2008. p. 323-327). Deterministic approaches like iterative and model-based method for segmenting large vessels in CTPA examinations have been proposed (see e.g. Ecabert O, Peters J, Walker MJ, Ivane T, Lorenz C, von Berg J, et al. Segmentation of the heart and great vessels in CT images using a modelbased adaptation framework. Medical image analysis. 2011;15(6):863-876). More recently, probabilistic approaches like deep learning based system used for the segmentation of large vessels (see e.g. Baskaran L, Al'Aref SJ, Maliakal G, Lee BC, Xu Z, Choi JW, et al. Automatic segmentation of multiple cardiovascular structures from cardiac computed tomography angiography images using deep learning. PloS one. 2020;15(5):e0232573). However, the previous studies have two main limitations; their datasets either (i) contain an insufficient number of CT examinations or (ii) do not reflect the characteristics of CTPA examinations obtained in daily routine. The deep learning-based approaches were trained and tested with a large number of CTPA examinations, but there is no information about dataset characteristics. To overcome these limitations, there is a need to train and test algorithms on large, unselected datasets containing representative artifacts and medical conditions often encountered in radiology practice such as pericardial fluid, pleural effusion, and pulmonary embolism. Accordingly is provided an automated measurement of morphometric and geometric parameters of mediastinal structures from CTPA examinations by a fully automated algorithm that could perform equally well or better than an experienced radiologist. Thus, an algorithm that automatically measures mediastinal cardiovascular parameters which can aid diagnosis of cardiovascular conditions such as PE and PH in nonelectrocardiogram (ECG)-gated CT pulmonary angiography (CTPA) examinations is provided.

The current disclosure provides solutions to the above-mentioned problems and drawbacks by providing an algorithm that automatically measures mediastinal cardiovascular parameters, aiding diagnosis of cardiovascular conditions such as pulmonary embolism and pulmonary hypertension in CTPA examinations, an automated method and system which could scan obtained medical images for abnormalities automatically. An automated computer aided method is thus provided, which detects these abnormalities automatically, and that indicates said results to the clinician. Accordingly, the present disclosure provides a fully automated CADe algorithm for measurement of morphometric and geometric parameters of large vessels in CT pulmonary angiography examinations, which may be used for detection and diagnosis of conditions or disorders in the subject being scanned. The present method includes an algorithm that is capable of performing the measurements automatically by automatic detection of automatic detection of the seed points, which allows for organ detection and extraction, which may be used to automatically identify abnormalities in the vessels, such as by determining the average size (diameter) of the ascending aorta and pulmonary trunk in obtained images from a subject automatically. A CT scan of the subject give rise to CT data, the CT data comprising information comprising an image stack forming a CT volume of images. The CT data may be analyzed using an algorithm implemented in a software tool, an algorithm for automatic detection of the seed points, which allows for organ detection and extraction using measurements of morphometric and geometric parameters of large vessels, such as the pulmonary trunk, ascending aorta and descending aorta. The method may for example provide an automatic method for determining a representative diameter, corresponding to the average diameter, of the pulmonary trunk or ascending aorta. The representative diameter may then be compared to one or more predetermined thresholds, and a diameter above one such predetermined threshold may be indicative of pulmonary hypertension or aneurysm. Thus, the result of importing the CT data into the computer software and running the algorithm of the software program/tool may then be an indication whether e.g. the size of the pulmonary trunk in the CT is above the threshold or not, i.e. indicate pulmonary hypertension or not. If the result indicates that pulmonary hypertension is present in the subject of the corresponding CT data, a clinician may further use the attained representative diameter as decision support for an eventual PE diagnosis, and/or the clinician may further evaluate the images of said subject based on the indication of pulmonary hypertension being present to perform an eventual diagnosis. The provided method of automatic detection of seed points, and detection and measurements of the large vessels, may be useful for other purposes as well, where such an automatic seed point detection and identification of large vessels may be useful. These uses are also covered by the present invention. Non limiting examples comprise surgical planning, calcification detection, aneurysm detection, dissections and stenosis detection, to name a few.

Medical imaging is the technique and process to create images of the patient's interior body for medical diagnosis and treatment purposes. Nowadays medical imaging systems are indispensable for modern medicine that is incorporated with radiology. One widely used such imaging technology is CT. CT technology is mainly based on X-ray physics, which uses combinations of scanned X-ray slices (corresponding to the virtual image in radiology) to produce axial (transverse) images. In view of the development and the technology of CT such as X-ray tubes and detectors, there are several systems in use: Standard (Conventional) CT, Helical (Spiral) CT (HCT), High-Resolution CT (HRCT), Multiple-row Detector CT (MDCT), Dual Source CT (DSCT) are some of the most known systems. Each of them is different from the other in terms of scanning time and scanning speed, image quality, and dose efficiency.

The working principles of CT scanners can be divided into four basic steps. Completing all these steps may be referred to as performing a CT scan, and the resulting data comprising an image stack forming a CT volume of images may be referred to as CT data. First, the patient lies on a motorized table. Second, while the table moves into gantry, X-ray sources rotate around the table. In the third step, detectors record patient's body snapshots, which are exited from X-ray beams. In the fourth and final step, all snapshots are sent to the computer to reconstruct into axial slices (Figure 1). In the final reconstruction step, each snapshot is turned into an image matrix (slice), which is made of numerical values. This matrix consists of elements, and these elements correspond to pixels. According to the slice thickness, each pixel has a volume and is referred to a volume element or voxel (Figure 2). The voxel size is calculated by:

VoxelSize = PixelSize x SliceThickness [Equation 1]

The matrix size and slice thickness change due to the CT scanner types and generations. In first generation CT scanners, the matrix size is equal to 80 x 80, and slice thickness is over 8 mm. Over the CT generations, matrix size has increased gradually (256 x 256, 512 x 512, and 1024 x 1024), and the slice thickness has decreased gradually to 0.5 mm. The Image quality and the image detail improves as the matrix size increases and the slice thickness decreases. Each pixel in the image matrix is a representation of X-ray density, which is evaluated from X-ray detectors. These densities are attenuation coefficients for the X-ray beam and called the CT number. To honor the inventor of CT, these CT numbers are also known as Hounsfield Units (HU), which are obtained by the attenuation coefficient. Each tissue in a CT scan has a different radiodensity, and it is calculated with regard to the radiodensity of distilled water and air. The radiodensity of distilled water at standard temperature and pressure (STP) is equal to 0 HU, and the radiodensity of air at STP is equal to - 1000 HU.

The Hounsfield Units calculation is given as:

HU = 1000 X - Equation 2] jiwater-jiair

For various reasons, some deformations and undesirable effects can occur in CT images. These distortions are known as CT artifacts. Due to the underlying reason of the artifact, CT artifacts can be divided into four categories: physics-based artifacts, patient-based artifacts, scanner (hardware)-based artifacts, and helical and multi-section artifacts.

Image quality (Q) in CTPA exams is affected by several aspects such as motion-breathing artifacts, streak artifacts, image noise, contrast concentration and lung parenchymal diseases. If each aspect can be scored on a point scale, then the image quality can be calculated as the sum of these scores (s) and formulated as follows: Q(s) = l i P^si [Equation 3] where ps represents the point scale of the aspect and N is the total number of aspects. However, each aspect has a different added factor (k) on image quality, yielding:

[Equation 4]

In the methods of the present disclosure, the aspects, point scales and added factors were determined subjectively by radiologist TF. Using equation 2 with the parameters of selected aspects with point scales and added factors and equation 3, the image quality score s can be determined as:

[Equation 5]

are added factors and mba, sa, d, Id, and n represent the aspects motionbreathing artifacts, streak artifacts, contrast concentration, lung disease, and image noise, respectively.

When performing medical imaging of the body, it is important to consider in which plane of the human body the scan is performed. The different planes of the human body may be divided into axial (transverse or horizontal), which divides the body into superior (upper) and inferior (lower) regions, sagittal (median) that divides the body into left and right regions and coronal, which divides the body into anterior (front) and posterior (back) regions. To create image series in CT, in the last reconstruction step of the CT scanning, the raw data, which is taken from a CT scanner, is sent to a computer to reconstruct into axial (transverse) slices. The data attained from said reconstruction may be referred to as "CT data" herein. From axial series, it is reformatted in multiple planes regarding slice thickness and planes. For the most part, radiologists look at axial, sagittal and coronal series in CTPA images at the same time to detect and diagnose PE.

The thorax is the major part of the body that located between the neck and the abdomen part of the body, and thorax anatomy in CT scans covers said region. The thoracic cavity contains the heart, lungs (left and right), trachea, esophagus and other organs, and it is supported by the thoracic wall. The trachea is a tubular organ that, besides the lungs, is an important member of the respiratory system. It is also known as windpipe. Its function in the respiratory system is to transport air to the lungs. At fifth thoracic vertebra (T5), the trachea is divided into two branches: left and right primary bronchus. This separation is also called as the bifurcation of the trachea. After the bifurcation point left and right primary bronchus enters the left and the right lung respectively. Likewise, primary bronchi are divided into secondary (lobar) bronchi, where bronchus enter the lungs. Afterwards, each secondary bronchus are divided into tertiary (segmental) bronchi. The mediastinum area is a space between the left and the right lung that divided into the superior and inferior part. Inferior part is subdivided into three parts: anterior, middle and posterior mediastinum. The Mediastinum contains heart, pericardium, pulmonary vessels (arteries and veins), trachea, thoracic aorta, esophagus and the other organs. The major pulmonary arteries inside of the mediastinum is the main and lobar pulmonary arteries, pulmonary trunk, and pulmonary valve. Located between the left and right lung is the heart which pumps blood into the circulation system through the vessels. The major vessels cover pulmonary arteries and veins, the superior and the inferior vena cava and the thoracic aorta (ascending aorta and descending aorta).

The pulmonary trunk or main pulmonary artery (PA) is a vessel that arises from the right ventricle of the heart, extends upward, and divides into the right and left pulmonary arteries that convey unaerated blood to the lungs. The pulmonary trunk is the solitary arterial output from the right ventricle, transporting deoxygenated blood to the lungs for oxygenation. The pulmonary trunk is approximately 50 mm long and 30 mm wide, and 29 mm width is often used as the cut-off of normal. It arises as a direct superior continuation of the right ventricular outflow tract, separated by the pulmonary valve. As it ascends it slants posteriorly and to the left of the ascending aorta. With the ascending aorta, the pulmonary trunk is invested in a common sheath of serous visceral pericardium, anterior to the transverse pericardial sinus. At the level of the transthoracic plane, the trunk emerges from the fibrous pericardium and divides into the longer right and shorter left pulmonary arteries in the concavity of the aortic arch, anterior to left main bronchus and to the left of the carina. The left coronary artery passes between the pulmonary trunk (on the left) and the auricle of the left atrium. The size of the pulmonary trunk is linked to pulmonary hypertension which may indicate PE, and may be determined using automated image processing.

Image processing is a subdiscipline of signal processing which deals with the processing of analog and digital signals. All kind of electronic snapshots are an example of digital images (digital signal) and is thus referred to as digital image processing. Simply, an image is a dimensional signal which can be represented as mathematical function that defined by f(x,y)=g for 2D and/f'x,/, z)=g for three-dimensional (3D) volume images, where x, y, and z are the coordinates of the smallest element of the image and g represents the gray level value of that element. In 2D images, this element is known as picture element (pixel). In 3D images, this element is referred to as volume element (voxel). In sum, pixel and voxel values are corresponded to gray level value of images at spatial coordinates.

Image processing and image analysis tasks are separated from each other concerning the outputs obtained. In an image processing task, the output is a different aspect of the input image. However, in an image analysis task, the output is valuable information that is extracted from the input image. Image acquisition is the first step of image processing where images are retrieved as a dataset from a source such as medical devices, satellites, cameras, as in the present disclosure where the dataset may consist of CTPA images. These images are retrieved from the CT scanner and formatted and recorded under DICOM standard, and this process is the image acquisition step of the image analysis task.

As shown in Figure 3 is an illustration of an image acquisition process in three steps, the image acquisition process resulting in CT data comprising an image stack forming a CT volume of images . In step 1, a CT scan of the subject is performed by the CT scanner (100) resulting in a number of slices or scans (15), and the subject's scans (15) are send to a computer (50) for reconstruction. In step 2, each scan (15) is turned into an image matrix, using a software in the computer (20) where a digital image is represented as a matrix function in a computer system that defined by 2551

[Equation 6]

120J for 2D images, where m and n are the sizes of the matrix, and each element in the matrix corresponds to the pixel value. In step 3, image matrices are reformatted and stored for example under DICOM standard on servers (25). The reconstructed, reformatted image matrices form CT data (20) stored on the servers (25), which may then be obtained for image analysis.

The retrieved images may undergo image enhancement. Image enhancement is the process of improving the appearance of images for further image processing and analysis steps or human viewers. Noise reduction, histogram equalization, contrast adjustment, image sharpening, logical operations and arithmetic operations are some methods of image enhancement that may be performed. Morphological image processing may be performed on binary images or grayscale images, which may be used to smooth contours of an object or to fill gaps.

The most crucial step in an image analysis problem is the segmentation task. Segmentation subdivides an image into meaningful regions to solve and analyze the problem efficiently. A typical segmentation task is carried out by several methods and approaches by using some characteristics of a region (or an object) such as color, density, pattern, and morphology. Different segmentation methods exist, such as thresholding mostly used for converting a grayscale image to a binary image, edge detection used for finding boundaries of an object by detecting these discontinuities, and region growing which has been used in the experimental part of this disclosure.

Region growing is a region-based image segmentation method widely used in medical image segmentation tasks. The first step in the region growing method is the selection of the seed points. To find adjacent points from these seed points, a region growing strategy is determined. This strategy could be such as connected neighborhood for binary images or pixel intensity for grayscale images. In the next steps, the region is grown iteratively depending on these strategies. The region growing method could be performed on 2D or 3D images. For 2D binary images, 4 or 8 connected neighborhood can be used. For 3D binary images, 6, 18, or 26 connected neighborhood can be used. In image processing, watershed refers to an advanced region-based image segmentation method that is inspired and improved based on geomorphology definitions. Watershed transform treats the grayscale image like a topographic map to find watershed lines and catchment basins. In the grayscale image, bright pixels represent high elevations, and dark pixels represent low elevations.

Feature recognition and classification is the final step of the all image processing and analysis problem. Detecting or recognizing an object is a feature recognition task, e.g. does this CTPA image show a pulmonary embolism? An object can be detected or recognized by some set of features such as edges, blobs, corners, and ridges. Dividing an object into the different group or classes is a task of classification, e.g. is this pulmonary embolism acute or chronic? Classification can be performed by three methods: statistical methods, structural methods, and hybrid methods.

Automatic interpretation of CT volume images to measure/determine for example a representative diameter of the pulmonary trunk or ascending aorta requires segmentation and extraction tasks such as airway tree extraction, vessel segmentation, and pulmonary trunk detection and extraction, after which lung CT volume images can be used for the computer-aided diagnosis (CAD) of various diseases such as pulmonary embolism detection.

The extraction of airway trees is one of the most difficult tasks in CT volume image segmentation tasks. The airway extraction task can be solved by several approaches, but the two main approaches are common; region growing based methods and mathematical morphology based methods. Vessel segmentation is one of the well-studied tasks in the CT volume images segmentation tasks. The vessel segmentation task can be solved by several algorithms such as pattern recognition techniques, model-based approaches, tracking-based approaches, artificial intelligence-based approaches, neural network-based approaches, and tube-like object detection approaches. However, three main approaches stand out; the hessian-based method, the region growing based method and the thresholding based method. Since the automatic interpretation of lung CTPA volume images requires several segmentation tasks, different methodologies are employed for each of the tasks.

Example operations

The proposed methods will now be described in more detail referring to Figure 4. It should be appreciated that Figure 4 comprise some operations and modules which are illustrated with a solid border and some operations and modules which are illustrated with a dashed border. The operations and modules which are illustrated with solid border are operations which are comprised in the broadest example embodiment. The operations and modules which are illustrated with dashed border are example embodiments which may be comprised in, or a part of, or are further embodiments which may be taken in addition to the operations and modules of the broader example embodiments. It should be appreciated that the operations do not need to be performed in order. Furthermore, it should be appreciated that not all of the operations need to be performed.

Figure 4 relates to a method, which may be performed in a computer or computing device (10). The method may be performed using a computer program/software, comprising computer program code which, when executed in the computer or computing device (10) causes the computer or computing device (10) to execute the methods of the present disclosure. Figure 4 illustrates a method, such as a computer-aided detection, CADe, method performed in a computing device, for locating large vessels in a subject undergoing a medical imaging procedure, such as a computed tomography (CT) scan. The located vessels may also be measured automatically, such as measuring the diameters of the large vessels for identifying possible abnormalities which may be indicative of a disorder or medical condition, such as measuring and determining a representative diameter of a pulmonary trunk and/or ascending aorta, wherein a dilated diameter may be indicative of PH and/or PE. The method and system may first automatically locate the three-dimensional positions of two anatomical landmarks, the carina of trachea and an apical level of the pulmonary valve (PV), to find seed points of vascular structures in the mediastinum (ascending aorta, descending aorta, pulmonary trunk). Then, the seed points for segmentation of mediastinal structures may be automatically placed using a heuristic approach, after which the structures may be segmented by applying image processing techniques, such as image enhancement, edge detection, gray scale segmentation, and 2D region growing methods. After segmentation, measurements may be done on the segmented regions.

An example embodiment of the invention comprises a first step of obtaining (S10) CT data from a consecutive CT scan of the thorax of the subject, the CT data comprising an image stack of slices forming a CT volume of images. Obtaining the CT data may be done by retrieving the CT data stored on a sever, retrieving, using the computing device, stored image matrices as CT data from a server to obtain the CT data. The CT data on the server may in turn be acquired by performing a CT scan of the subject using a CT scanner (100) to retrieve a number of slices, sending the retrieved slices to a computer/computing device for reconstruction, turning each slice into an image matrix using a software in the computing device, and reformatting and storing the image matrices under DICOM standards, or other applicable standard, in servers by sending the DICOM folder to e.g. PACS servers for storing the patient data, from which server(s) the CT data may then be retrieved, by e.g. other computing devices. The image stack of consecutive CT scans may be computed tomography pulmonary angiography (CTPA) volume images.

The next step may comprise preparing (Sil) the dataset for segmentation tasks. The CT dataset may be prepared for segmentation tasks by first calculating (Sila) a linear scale value. Stored values (SV) in DICOM files are attenuation coefficients. Therefore, these stored values must first be converted to a linear scale (Hounsfield Unit scale for CT volume images) by the following equation:

HU = m x SV + b [Equation 7] where m is the slope, SV the pixel value and b the intercept. Thus, in a first step, every voxel in the CT Data (original linear attenuation coefficient) is converted to the Hounsfield Units scale using the equation 7, where the slope and intercept values are found in the DICOM header file tag (0028, 1053) and (0028, 1052), respectively.

When operating, the developed system of the present disclosure first locates the three dimensional positions of two anatomical landmarks, carina of trachea and an apical level of the pulmonary valve, to find seed points of mediastinum structures (ascending and descending aorta, pulmonary trunk). These anatomical landmarks are detected automatically by the system scanning CTPA volume images in the cranial to caudal direction. Hence, knowing the scanning direction is crucial to detect the anatomical landmarks. Further, the orientation of the subject in the CTPA image with respect to the x-axis of the image plane is a key factor to automatically and accurately detect the carina of trachea, the apical level of the pulmonary valve, the ascending and descending aorta and the pulmonary trunk.

Thus, the next step is obtaining (Sllb) information indicating a scanning direction of the CT scan. In a CT scan procedure, the subject can be scanned in two different directions; in cranial to caudal (head to tail or superior to inferior) direction or caudal to cranial (tail to head or inferior to superior) direction depending on the CT examination procedure. If the scanning direction is cranial to caudal, then the first images in the 3D CT image data are the superior part of the thorax. If the scanning direction is caudal to cranial, then the first images in the 3D CT image data are the inferior part of the thorax. The scanning direction can be calculated by comparing the first and last CT slice position (patient) information in the DICOM header file tag no (0020, 0032). For example: The z-axis value of the first slice according to DICOM header file tag (0020, 0032) is -280.75 The z-axis value of the last slice according to DICOM header file tag (0020, 0032) is 48.75 The smaller z-axis value of the first and last slice indicates the caudal part of the patient. In this example, the caudal part starts with the first CT slice and the scanning direction is caudal to cranial.

After the scanning direction has been obtained, determining (Sllc) an orientation of the subject in the scanner is performed, wherein determining said orientation of the subject comprises calculating the subject's orientation in an image of the image stack as an angle a (also referred to as the "a -value") between a major axis of the image and an x-axis of an image plane; and on a condition that the angle a is larger than 17 °, performing an additional calculation of the subject's orientation based on caudal slices of the CT scan.

The orientation of the subject is the angle between the major axis of the subject's image and the x-axis of the image plane (Figure 5). Figure 5 illustrates orientation of a CT slice. The upper lines of a) b) and c) represent the major axis of the patient and the lower lines represents the major axis (x-axis) of the image, (a) Both lines are parallel to each other. Therefore, the angle between the upper and lower lines (the orientation of the CT slice) is near 0, which is required for proper functioning of the automated system. In (b) the angle between the lines is 21 degrees (alpha = 21 degrees). In order to detect anatomical landmarks (carina of trachea, level of the pulmonary valve) and major mediastinum structures (ascending and descending aorta, pulmonary trunk) precisely, the major axis of the patient must be rotated 21 degrees respect to the major axis of image. The rotated CT scan is represented in (c).

To achieve accurate detection, the major axis of the subject in the image must be parallel to the x-axis of the image plane. The first image in a CT volume may contain information about the volume rather than images of the subject. Therefore, the cranial orientation of the subject is first computed using the second image of the CT volume images (Figure 6 and Figure 7). Figure 6 show algorithmic steps of a first patient orientation calculation method. First, acquire the second image from the CTPA volume images, then perform image thresholding, the determine the largest component in the image, and last, calculate patient orientation. Figure 7 show algorithmic steps of a second (alternative) patient orientation calculation method. Figure 7 illustrates a method developed by the present inventors to calculate patient orientation more precisely based on the caudal slices of the CT exam, as their orientation was found to correlate closer to the orientation of the carina of trachea and pulmonary valve. In a first step of Figure 7A, the second to last image from the CTPA volume images was acquired. To know the patient orientation in the caudal part, it is sufficient to examine a single image from the caudal part. It was observed that doing the calculation with anyone of the last slices in the CTPA volume is adequate for determining orientation correctly. Therefore, the second to the last slice from the volume was selected. After acquiring the image, image thresholding was performed as a second step of Figure 7A. The image were thresholded over -300* HU in order to detect the largest component. The threshold of -300HU was chosen based on literature and empirical observation. Then the largest component in the image was determined (Figure 7A step 3). Using 2D connected component analysis (MATLAB R2019b bwconncomp function), one may determine the largest component. Next, perform a morphological opening operation to clean the boundaries of the component. In step 4 of Figure 7A, find region boundaries of largest component. The region boundaries are obtained by applying a boundary mask function (MATLAB R2019b boundary mask function). In step 5 of Figure 7A, divide largest component into four regions. Divide the largest component into 4 regions by drawing perpendicular dividing lines starting from the center of the image. In step 6 of Figure 7A, extract the 3rd region. The method continues in Figure 7B, with step 7 of extracting the 4th region. Empirically, the patient orientation can be computed from information in regions 3 and 4. Hence, these regions are extracted from the image. In step 8, extract 20 pixel area from the bottom upwards along the y axis. To decide which one of region 3 or 4 that provides more information about the orientation, segment the region bounded by 20 pixels from bottom to top on the y axis and all pixels on the x axis. Next, count the number of white pixels in this segmented region. In step 9, determine direction of tilting. Let the sum of all white pixels for segmented area in 3^rd and 4^th region be s3 and s4, respectively. Assume that the scalar quantity of pixels in this region indicates which direction the object is tilted. Then, the CT scan is tilted to the right (region3) if s3 > s4. If s4 > s3, the CT scan is tilted to the left (region4). In step 10, apply polynomial curve fitting to the 3rd region, and in step 11, apply polynomial curve fitting to the 4th region. In last step 12, choose tilting direction and angle, and apply rotation to the image. If we find the straight lines which are the respective best fits to the curve in region 3 and 4, we can calculate the angle of these lines with respect to the x axis. The best fitted line for a given function (a series of datapoints) can be calculated by polynomial curve fitting. After calculating angle alpha, the entire CT volume is rotated by alpha.

Subject orientation may vary along the superior to inferior direction in the CTPA study because of scoliosis, subject movement during the examination, the subject lying down obliquely on the CT table during the examination, or other reasons. Curvature on just few slices on superior or on inferior part of subject's image does not interfere with the proper functioning of the developed algorithm. However, if the angle the a is larger than 17 ° in any of the first slices (top 10 slices) in the CTPA study, such as the second slice, an additional orientation test is performed on the caudal part of the patient. When computing patient orientation with respect to the x-axis, the anatomy of the patient may preclude correct orientation (Figure 8). Figure 8 shows obtaining accurate orientation of the CT scan. The orientation of the original image in (a) is 27°. Thus, the image should be rotated 27°. If the orientation is calculated with the method proposed in step 3, the result is 12° and the rotated CT scan can be seen in (b). The incomplete rotation interferes with automatic detection of carina of trachea and the pulmonary trunk. However, with the method proposed in the method above, the calculated orientation is 30° (c) which is sufficiently accurate for proper functioning of the developed system.

Therefore, a new simplified method to calculate patient orientation more precisely based on the caudal slices of the CT exam was developed by the present inventors, as their orientation was found to correlate closer to the orientation of the carina of trachea and pulmonary valve. One way to find the orientation of the patient is to calculate the orientation of the spinal canal. To do this, the spinal canal must first be detected in the CT volume, and the entire spinal canal must be segmented correctly. However, calculating the orientation of the patient in this way is computationally challenging and hard to validate. However, with the method proposed one can easily calculate the orientation. The method comprises acquiring the second image from the CTPA volume, perform image thresholding, determining the largest component in the image and calculating the subject orientation. Instead of calculating the CT scan orientation on every image, it is sufficient to calculate it only in a single image in the first place. It was observed that basing the calculation on any one of the first slices in the CTPA volume gives the best results, such as the second slice. In an example, the image is thresholded over-300 HU in order to detect the largest component (surface area) in the next step where the largest component is determined using 2D connected component analysis (e.g. MATLAB R2019b bwconncomp function). The image is input into the region props function (MatlabR2019bimageprocessingtoolbox) which calculates the a value. The a cutoff value 17 (the angle a is larger than 17 °) was empirically determined.

Thus, if the new simplified method yield an angle a which is equal to or less than 17°, i.e. not larger than 17 °, it is determined that the attained patient orientation is sufficiently accurate and the CT scans of the CTPA exam are rotated according to the angle a. However, if the new simplified method yield an angle a which is larger than 17°, it is determined that the attained angle a may be insufficiently accurate, and thus an additional orientation test according to the more complex method is performed, an additional calculation of the subject's orientation based on caudal slices of the CT scan to attain an updated angle a.

Thus a step of performing (Slid) an additional calculation of the subject's orientation based on caudal slices of the CT scan to attain an updated angle a is performed. The additional calculation comprises acquiring the second to last image from the CTPA volume images, performing image thresholding, determining the largest component in the image, finding region boundaries of the largest component, dividing the largest component into four regions, extracting the third region, extracting the fourth region, extracting a 20 pixel area from the bottom upwards along the y-axis, determining direction of tilting, applying polynomial curve fitting to the third region, applying polynomial curve fitting to the fourth region and choosing tilting direction and angle, where the angle a is determined from the angles calculated by the polynomial curve fittings of region 3 and 4. The attained angle a from the more complex method, referred to as an updated angle a, is then used for applying rotation to the CT volume.

In summary, to decide how to rotate a CT scan or CT volume of images (CT data) for image processing and analysis, the angle a in the cranial part of the subject is first checked using the simplified method above, and if the cranial angle a > 17 °, an updated angle a is computed in the caudal parts of the patient by the more complex method above. Finally, all CT scans of the CTPA exam are rotated according to the angle a or updated angle a. Thus, after the determining the orientation of the CT scan and obtaining the angle a; a step of correcting (Slle) scan curvature by rotating the CT scan according to the computed angle a, or if present, the updated angle a, is performed. The CT data has now been prepared for segmentation tasks.

After these steps preparation steps, the image volume is ready for the computer aided method of locating and measuring large vessels. The first step of said method relates to locating (S12) carina trachea in the CT volume, which comprises locating (S12a) the trachea of the subject in the CT volume of images, detecting (S12b) if any tracheal intubation is present, segmenting (S12c) the airways of the trachea region; and identifying (S12d) the optimal carina location for the trachea at carina level as the carina trachea. Locating (S12a) the trachea in the CT volume comprises: downscaling the CT volume to half its size, determining a first volume of interest (VOI1), segmenting air areas in VOI1, generating a pool of potential trachea candidates (TCP1) in VOI1, and on a condition that there is one 3-dimensional, 3D, component in the TCP1: locating the trachea as the 3D component, and on a condition that there is more than one 3D component in the TCP1: determining a second volume of interest (VOI2), segmenting air areas in VOI2; generating a second pool of potential trachea candidates (TCP2) in VOI2, combining TCP1 and TCP2 to find common 3D components, and locating the trachea as the 3D component spanning the largest total number of CT images of the combined TCP1 and TCP2. When the trachea is located, it is possible to select a slice from the CT volume of images where the trachea is located as a seed point S_t for a subsequent extraction of the trachea.

As an example, in a majority of CT exams, the trachea can be located as a tubular air-filled structure of lower HU density than surrounding lung tissue. However, in a minority of CT exams, there are artefacts in the CT planes such as beam hardening, craniofacial structures or upper airway structures that resemble the trachea. Finding the trachea in presence of such artefacts have received little attention in the prior art. To accelerate computation, the CT volume I (x, y, z) is downscaled to half its size (S12a, step 1) with respect to the x- and y- axes but not the z-axis: z) . In some CT exams, the trachea adjoins the right lung

before the carina of the trachea because of the curvature of the trachea. This prevents localization and segmentation of the trachea.

In an example, in an empirical observation, such adjoining did not occur in the cranial 15% of CT slices in 271 CT exams. To correctly localize the trachea in 3D, the system therefore generates two consecutive volumes of interest where trachea candidates are generated and assessed. First, the cranial of the two consecutive volumes of interest, volume of interest 1 (VOI1) is determined by extracting the CT slices starting from the second slice (S2) of the CT exam up to s_n in the cranial to caudal direction, where s_n = 15 % of the total number of slices in the CT scan (S12a, step 2). Air areas (lungs, airways, and artefacts that resemble air) are segmented in VOI1 by applying the following steps (S12a, step 3) to each slice of VOI1: The slice is thresholded over -300 HU, and 2D connected component analysis is applied. The component with the largest area is designated as the thoracic cavity. A morphological floodfill operation is applied. To get air filled areas, a logical AND operation is applied the output of step (S12a, step 1) and the output of applying the morphological flood-fill operation.

Next, a pool of potential trachea candidates (TCP1) in VOI1 is generated (S12a, step 4). The average HU density of VOI1 is calculated. All slices in VOI1 are then thresholded with this average density. Empirically, the area of the trachea was found to be within the range 10- 1200 pixels in the CT slices of 271 CTPA exams. Therefore, components with an area < 10 or > 1200 pixels are excluded in each slice. To remove false positive trachea candidates, volumes < V voxels are excluded where V is calculated by:

V = ( | s_n — 1 +1) *7 [Equation s] where s_n and S2are calculated in above is an empirically determined coefficient to calculate likely minimum volume of trachea candidates.

When the above steps (S12a, steps 1-4) are applied to the CT volume, most often more than three 3D components (trachea, esophagus, right lung, and left lung) are obtained, but in some cases only two 3D components (trachea, right or left lung) are obtained. If there is only one 3D component in the TCP1 as a result of above steps, this component is designated as the trachea. However, if there is more than one component in the TCP1, it is not obvious which component is the trachea, the right lung, the left lung or the esophagus. Therefore, to recognize the trachea correctly and automatically the following steps are applied:

S12a, step 8: A second volume of interest (VOI2) to search for the trachea is defined. The volume starting from the last slice of VOI1 (s_n) up to slice st which is calculated by

Sfc =m-m*0.25 [Equation 9] where m is half of the total number of slices in the CT volume and the empirically determined coefficient 0.25 is used to calculate the total number of slices to be included in VOI2 based on empirical observation from 271 CT exams. The VOI2 will contain the right and the left lung, except when one of the lungs is completely collapsed. In step S12a, step 9, segmentation of air areas in VOI2 is performed. Air areas (lungs, airways, and artefacts that look like air) are segmented in VOI2 by applying the following steps to each slice of VOI2:

1. The slice is thresholded over -300 HU.

2. 2D connected component analysis is applied.

3. The component with the largest area is designated as the thoracic cavity.

4. A morphological flood-fill operation is applied.

5. A logical AND operation is applied the output of step 1 and the output of step 4 in order to get the air filled areas.

After, step S12a, step 10, generating a second pool of potential trachea candidates, TCP2, in VOI2 is performed. A pool of potential trachea candidates (TCP2: trachea candidates pool 2) is generated for VOI2 by the following steps:

1. The average HU density of VOI2 is calculated.

2. The average HU density computed in 1 is multiplied by 0.05 which is a HU density enhancement factor determined by empirical observation.

3. All slices in VOI2 are thresholded with the HU density computed in 2.

4. Areas < 10 and > 1200 pixels in each slice are excluded.

5. Volumes < V voxels are excluded, and V is calculated by follow:

V =( | Sk-Sn | +1) *7 [Equation 10] where Sk is calculated in S12a, step 8, and s_n is calculated in step S12a, step 3 and the constant 7 is an empirically determined coefficient to calculate likely minimum volume of trachea candidates.

The next steps, S12a steps 11, 12 and 13, combining TCP1 and TCP2 to find common 3D components, locating the trachea as the 3D component spanning the largest total number of CT images of the combined TCP1 and TCP2, and selecting a slice from the CT volume of images where the trachea is located as a seed point St for a subsequent extraction of the trachea is performed. Choosing the right component that is the trachea comprises: To mark some 3D components to potentially be a trachea (TCP1) one may limit the search space (VOI 1) in S12a steps 1-4. However, it may not be known which component in the TCP1 is the trachea. Therefore, in steps S12a, step 8 to S12a, step 10 a second search space is generated (VOI2) from the CT scan, which is consecutive of VOI1 in cranial to caudal direction, to limit the trachea candidate pool. Anatomically in VOI1+VOI2, the longest 3D component in cranial to caudal direction is the trachea. In the following steps, TCP1 and TCP2 are combined to detect the longest component:

1. Each 3D connected component in the TCP1 is labeled using the bwlabeln function in MATLAB R2019b.

2. Each 3D connected component in the TCP2 is labeled using the bwlabeln function in MATLAB R2019b.

3. To find common components between TCP1 and TCP2, a logical AND operations is applied to the outputs of step 1 and step 2 at slice , which is the only common slice between the two volumes.

4. To obtain all components related to TCP2, a logical AND operation is applied the output of step 1 and step 3.

5. To obtain all components related to TCP1, a logical AND operation is applied the output of step 2 and step 3.

6. The outputs of step 4 and 5 are combined by arithmetic addition.

7. The longest component, defined as the component spanning the largest total number of CT slices in the combined volume of 6, is designated as the trachea.

8. For further steps, one must choose one of the CT slices along the z axis where the trachea is located by : s_t =round( min (p) + max (p) /2) [Equation 11] where is a matrix containing positions of the trachea along the z-axis.

Next step after locating the trachea is detecting (S12b) if any tracheal intubation is present. Some patients undergo CT examination with tracheal intubation. The insertion of an endotracheal tube creates a different appearance of the trachea in CT examinations and interferes with automatic detection of the carina of trachea. To overcome this problem, tracheal intubation needs to be detected by the system. If the trachea was successfully located in the previous step, tracheal intubation may be detected by the following steps; S12b step 1. The center point of the trachea is calculated by the bounding box or minimum bounding rectangle method. As the trachea appears as a circular object in a transverse CT slice, the center point of the smallest rectangle containing the trachea is the same as the center point of the trachea.

S12b step 2. Since the HU density of the inserted tube is high, the CT image at slice t from S12a steps 11, 12 and 13, is thresholded over 300 HU.

S12b step 3. Areas < 15 and > 100 pixels in the binary image of step 2 are excluded.

S12b step 4. 2D connected component analysis is applied to the binary image of step 3.

S12b step 5. The center points of all distinct 2D connected components in the binary image are calculated. If a tracheal intubation exists, then the endotracheal tube is represented by one of these 2D connected components.

S12b step 6. The endotracheal tube is located inside of the trachea. Therefore, the center points of the tube and the trachea must be close to each other. To find the closest components to the trachea, the Euclidean distance of component center points to the center point of the trachea is calculated:

[Equation 12] where n is the total number of components and x and y are the center points of the trachea and the components respectively.

S12b step 7. If the distance of the center points of 2D connected components to the center point of the trachea is < 11 pixels (corresponding to approximately 5,5 mm²) then these components are marked as potential endotracheal tubes. In some cases, calcification can occur around the trachea and the HU density and morphology of the calcification may appear similar to an endotracheal tube. Therefore, additional filters are required to ensure that the detected component is an endotracheal tube. The average HU density of the 2D component is calculated, and if >900 HU the component is designated as an endotracheal tube. While the HU density of a calcification may be >900 HU in theory, any calcification around the trachea have not been noticed where the HU density of the calcification is > 900 HU in areas > 15 pixels (corresponding to approximately 7,56 mm²) in a large dataset of 271 CTPA examinations. Once an endotracheal tube is detected inside the trachea, one need to find the CT slice where the intubation does not exist in order to track tracheal areas up to the carina of trachea accurately. The location of the new CT slice must be between the end of tracheal intubation and the carina trachea. One may find this location by the following steps;

S12b step 8. The 3D region growing method is applied to the CT scan in order to extract the tracheal intubation. The seed points for the 3D region growing method are generated from the step 5.

S12b step 9. The last CT slice number of the tracheal intubation in the cranial to caudal direction is obtained. This location set as the preferred CT slice number ( I) for tracking the trachea up to the carina trachea.

The next step in detecting the carina trachea, after locating the trachea and identifying any eventual tracheal intubation and to compensate for said eventual tracheal intubation, is segmenting (S12c) the airways of the trachea region, to be able to identify the optimal carina location.

Since the trachea is located before the carina in cranial to caudal direction in the CT stack I(x_t, yt, z_t) , it is possible to track the trachea up to the bifurcation point by comparing trachea regions slice by slice. To continue the tracking process, one need to segment the trachea region. One can segment the first trachea region by the 2D region growing method. The inputs of the region growing method are the binary image and seed points. First, the CT slice I zt is thresholded over -700 HU to acquire a binary image (S12c step 1). With this thresholding process one may only focus on airways voxels. While a threshold of -300 HU can be used instead of -700 HU, a threshold of -700 prevents the lung and airway areas from joining each other. Second, the location of the trachea in the 3D stack I t yt zt obtained from step 5 (if intubation has been detected, the location of the trachea was obtained from previous steps) will be designated as seed points (S12c step 2). One may then apply 2D region growing using these seed points ( t and yt) to segment the trachea region.

Prominent morphological changes in the trachea regions between two consecutive CT slices indicates a bifurcation point. Thus, once the trachea has been segmented, one may move to next CT slice I zt in cranial to caudal direction. One may threshold the new CT slice over - 700 HU (S12c step 3) and a logical AND operation is applied to this new binary image with the segmented trachea region (S12c step 4). The 2D region growing method is applied to the output of previous step (S12c step 5). The area of the segmented trachea is calculated as the total number of pixels in the segmented 2D component (S12c step 6). Empirically, the area of the airways in 2D was smaller than 1750 pixels before the bifurcation of the trachea in 271 CTPA exams. As thresholding at -700 HU is insufficient for distinguishing the trachea from the lung area, and basic image segmentation methods such as image thresholding are incapable of segmenting the trachea and lungs areas separately. If the calculated area exceeds 1750 pixels, the trachea is adjacent to the lungs. Therefore, one may apply watershed transform to segment the trachea from the lungs (S12c step 7). If the calculated area is smaller than 1750 pixels or once watershed transform is applied to the image, the trachea region in the previous CT slice I zt is subtracted from the newly segmented trachea region in the current CT slice I : : zt (S12c step 8).

So, a new component has been obtained. The difference between the newly segmented trachea region and the trachea region in the previous CT slice is the right main bronchus or the left main bronchus. To be sure about that the component, which is the difference of two images, is the one of the main bronchi, one may look for two prominent morphological changes. One is the mean densities of the components (component 1 and 2 represent the segmented trachea region in the current CT slice and the difference of two images, respectively) and the other is the size of the components. First, the mean densities of the components are calculated (S12c step 9). If the absolute value difference of the average density values of the components is less than 400 HU one may look for below condition in order to be sure that the difference of two images is the one of the main bronchi:

[Equation 13] where, Ac, At is the area of the component 2 and the area of the trachea region in previous CT slice respectively. And the k is the divide factor which is calculated as follow: }

[Equation 14]

The cut-off value 400 was determined by empirical observation in 271 CT exams. Once the main bifurcation on the trachea is detected, two components are obtained. The one is the right main bronchus and the other is the left main bronchus. One may then calculate the center points of these components and designate the left and the right main bronchi by comparing y-axes of the mass center of the components.

In a next step, S12d, identifying the optimal carina location as the carina trachea is performed. The step comprising: identifying the optimal carina location by tracking the trachea slice by slice from the bifurcation point in cranial to caudal direction, calculating the distance between left main bronchus and right main bronchus for each respective slice until reaching a slice, Soc, where the calculated distance is equal to or greater than 0,75 cm, and selecting said slice Socas the optimal carina location; and identifying the selected optimal carina location as the carina trachea. First one tracks up to the bifurcation point, then track again until the measured difference is reached.

Once the left and the right main bronchus is detected at the bifurcation point, one may track and segment the left and the right main bronchus, slice by slice, in cranial to caudal direction. In every slice, the left extrema points of the left main bronchus and the right extrema points of the right main bronchus are calculated. Subtract these points from each other, and if the result of the subtraction is equal or larger than 0.75 cm then designate this slice as an optimal carina point otherwise continue tracking the left and the right main bronchus slice by slice in cranial to caudal direction until the optimal carina point is reached.

The carina is a ridge of cartilage in the trachea that occurs between the division of the two main bronchi, thus the carina level consists of several CT slices. The "optimal carina location at carina level" is a representation of said carina level as a single CT slice, which is identified as the carina trachea in the method above. Since the carina trachea is a 3D structure, one need to designate a slice which represents the carina trachea in 2D. By empirical observation, the CT slice where the distance between the left main bronchus and the right main bronchus is at least 0.75 cm was designated as the optimal carina location I : : zc . This is performed by taking differences of vertical coordinates (y-axes) of the left most pixel of the right main bronchus (Pr) and the right most pixel of the left main bronchus (Py).This difference calculation operation is done in each slice starting from the where bifurcation occurred until the | Pr- Py | is equal or greater than 0.75 cm.

After step S12 of locating carina trachea in the CT volume, step S13 of detecting descending aorta at carina level is performed comprising: creating a first search space originating from an artificial line that cuts the x-axis of the slice at an angle of 22,5 -75 ° (such as 30 °, 40°, 50° or 60 °) downwards clockwise from the center of left main bronchus comprising lines having an angle larger than 2,5° and smaller than 15° between each other forming rays from center of origin, or an angle larger than 5° and smaller than 10° between each other, or an angle larger than 7° and smaller than 8° between each other, or an angle of 7,5° between each other. Creating a second search space originating from an artificial line that cuts the x-axis of the slice at an angle of 60 -120 ° (such as 60 °, 80 °, 100° or 120 °) downwards clockwise from the center of the left main bronchus comprising lines having an angle larger than 2,5° and smaller than 15° between each other forming rays from center of origin, or an angle larger than 5° and smaller than 10° between each other, or an angle larger than 7° and smaller than 8° between each other, or an angle of 7,5° between each other. Identifying if any one of the rays pass through the descending aorta by going ray by ray and: collecting all pixels on the ray, thresholding said pixels over - 1 HU to acquire binary image, applying 2D connected component analysis to identify largest connected pixels on the ray, calculating the total number of pixels on the segmented component, and if the total number of pixels is greater than 15, calculating mean x and y coordinates of those pixels, separately, in order to reduce to one data point, and identifying said data point as the descending aorta at carina level, PdesAo.

The descending aorta is located caudally to the left main bronchus. Since the left main bronchus was detected in a previous step, one may start from a line drawn from the center of the left main bronchus and that cuts the x-axis at an angle of between 22.5-75 degrees downwards clockwise. Thereafter, one may draw 8 lines, 100 pixels long, with an angle of between 2,5-15 degrees between each other (preferably around 7.5 degrees between each other) (Figure 9) to find the descending aorta (SSI: search space 1 Figure 9, left). However, in some cases, the descending aorta is shifted toward the right side. Therefore, one needs to shift the search space to the right. One may create a second search space (SS2: search space 2, Figure 9 right) by starting from an artificial line that cuts the x-axis at an angle of 60-120 degrees (preferably around 60) downwards clockwise from the center of the left main bronchus, and drawing 9 lines, 100 pixels long, with an angle of preferably 7.5 degrees between each other.

Once one has generated artificial rays, one may apply following steps ray by ray in order to locate the descending aorta. Those generated rays are used for taking pixels from tissues it passes over. One may apply the following procedures to see if the rays pass through an object like as the descending aorta. First, one may collect all pixels on the ray. Second, one may threshold those pixels over - 1 HU to acquire binary image. Third, one may apply 2D connected component analysis to find largest connected pixels on the ray. Fourth, one may calculate the total number of pixels on the segmented component. If the total number of pixels is greater than 15, then one may calculate mean x and y coordinates of those pixels, separately, in order to reduce to one data point. If the above criteria are met only by one ray in other word if one only has one x and y coordinates (PdesAo), one may assume that the pixel in this coordinate belongs to the descending aorta. Otherwise, If the above criteria are met by more than one ray then apply k-means clustering algorithms to the only x coordinates (since the rays radiate on x-axis) in order to cluster rays into two groups. After clustering the x coordinates, one may compare the number of elements between the clusters. The cluster with the most elements is selected. Then the x coordinates in this cluster and the y coordinates of these x coordinates are averaged to reduce one data point (PdesAo) . Once one have generated artificial rays, one may apply following steps ray by ray in order to locate the descending aorta. Those generated rays are used for taking pixels from tissues it passes over. One may apply the following procedures to see if the rays pass through an object like as the descending aorta. First, collect all pixels on the ray. Second, threshold those pixels over - 1 HU to acquire binary image. Third, apply 2D connected component analysis to find largest connected pixels on the ray. Fourth, calculate the total number of pixels on the segmented component. If the total number of pixels is greater than 15, then calculate mean x and y coordinates of those pixels, separately, in order to reduce to one data point. If the above criteria are met only by one ray in other word if there is only one x and y coordinates PdesAo one may assume that the pixel in this coordinate belongs to the descending aorta. Otherwise, If the above criteria are met by more than one ray then one may apply k-means clustering algorithms to the only x coordinates (since the rays radiate on x-axis) in order to cluster rays into two groups. After clustering the x coordinates, compare the number of elements between the clusters. The cluster with the most elements is selected. Then the x coordinates in this cluster and the y coordinates of these x coordinates are averaged to reduce one data point (PdesAo) . Once designated a reference point (PdesAo) which is possibly belongs to the descending aorta, one may apply the following methods of Gaussian filtering, Eigen values of Hessian matrix, Canny edge detection, gray scale segmentation, morphological operations, 2D region growing, morphological operations, in order, to segment and track the descending aorta slice by slice. Once obtaining the descending aorta, one may draw a 1 cm circle, the center point of which is the same as the mass center of the descending aorta. Then calculate the average HU density of this circle.

The next step, S14 detecting the ascending aorta, is performed by: detecting (S14a) aortic arch by tracking the descending aorta up to the aortic arch in caudal to cranial direction, detecting (S14b) the ascending aorta by tracking the upper part of the aortic arch to the ascending aorta in cranial to caudal direction to locate a first slice of the ascending aorta.

In the previous step, the descending aorta at the carina level was located /(:, : ,z_c).One can then track the descending aorta up to the aortic arch /(:, : ,z_aa) in the caudal to cranial direction. Prominent morphological changes in descending aorta regions between two consecutive CT slices indicates that one have reached the aortic arch. Therefore, one needs to segment all descending aorta regions /(:, : ,z_c: z_aa) between the carina of the trachea and the first prominent appearance of aortic arch in the CT stack. Segment and compare the descending aorta regions slice by slice by using the following methods: Gray scale segmentation, Anisotropic diffusion filtering, Eigen values of Hessian matrix, Canny edge detection, Morphological operations, and 2D region growing.

Figure 10 shows the aortic arch, where the area delineated by an oval line is the aortic arch in the CT slice. The vertical line is an artificial line, which divides the aortic arch into two parts; the upper part (indicated by dark triangle) and the lower part (indicated by light triangle). Once the aortic arch has been located, one can track the upper part (Figure 10) of the aortic arch in the cranial to caudal direction in order to reach the ascending aorta. Anatomically, tracking the upper part of the aortic arch for a few CT slices (between 2 to 7 slices) is sufficient for reaching the ascending aorta. Thus, tracking the upper part of the aortic arch for 3 slices starting from the first descending aorta regions in the CT stack in cranial to caudal direction. The upper part of the aortic arch is segmented by taking the difference of the first appearance of aortic arch in the CT stack from the first descending aorta regions in CT stack which are extracted in the previous steps. One may then track the upper regions of the aortic arch slice by slice in cranial to caudal direction by using the following methods: Gray scale segmentation, Anisotropic diffusion filtering, Eigen values of Hessian matrix, Canny edge detection, Morphological operations, and 2D region growing.

When the ascending aorta has been located, the program may perform a next step of segmenting and measuring the ascending aorta. The measurements may also be performed in a later stage after all the major vessels (ascending aorta, descending aorta and pulmonary trunk) have been located/detected, in which the method continues with detecting the pulmonary trunk in a next step.

Thus, optionally, in a next step, measurements for determining (S15) a representative diameter of the ascending aorta is performed, comprising: segmenting (S15a) the ascending aorta, tracking (S15b) the ascending aorta slice by slice from the first located slice to the carina trachea, wherein a diameter of the segmented ascending aorta is measured in every slice, and calculating (S15c) a representative diameter of the ascending aorta as a mean of the measured diameters from each slice.

In the previous step, the ascending aorta in the CT stack I (: ,: ,z_aa+2) was located. The next step is to segment all the ascending aorta regions between the first located slice I (: ,: ,z_aa+2) and the carina trachea I (: ,: ,z_c). Within this volume of interest / (: ,: , z_c'.z_aa+2)., one can segment and track the ascending aorta slice by slice using the following methods: Gray scale segmentation, Anisotropic diffusion filtering, Eigen values of Hessian matrix Canny edge detection, Morphological operations, and 2D region growing. In every slice, one may measure the diameter of the segmented ascending aorta. Finally, calculate the diameter of ascending aorta as the mean of these diameters. The diameter of the ascending aorta may be measured and calculated by using 'EquivDiameter' property of the 'regionprops' function of MATLAB R2019b Image Processing Toolbox. The returned scalar of this function is multiplied by pixel spacing value and the result is potential diameter. The representative diameter is calculated as a mean (+ standard deviation) from the measured diameters of each slice.

In a next step, detection (S16) of the pulmonary trunk is performed, comprising creating (S16a) a rectangular search space adjacent to the left lateral side of the ascending aorta location in the CT image stack covering the candidate of the pulmonary trunk, the pixel location of the rectangular being calculated as ((eAot - 15) :eAob , (eAol +5) : (eAol +70)), where eAot is the top most and eAob the bottom most pixel of the ascending aorta region; and detecting (S16b) the pulmonary trunk within said rectangular search space.

After determining the rectangular search area, calculate the mean pixel values (HU density) of all pixels with a pixel value greater than 45 HU. If the mean pixel values smaller than 125 HU then we threshold the rectangular search area over 45 HU otherwise we threshold the rectangular search area over the values subtracted 50 HU from the mean pixel value (mean HU - 50HU). After the thresholding process, mask the thresholded image with the original image. Further, one may enhance the rectangular search area by using the following methods: Anisotropic diffusion filtering, eigenvalues of Hessian matrix, Canny edge detection, morphological operations, and 2D connected component analysis. After the above procedures, if the calculated area is smaller than or equal to 350 pixels, repeat the above procedure again on the next slice in cranial to caudal direction. Otherwise, if the calculated area is greater than 350 pixels, designate this component as a potential part of the pulmonary artery system.

If the marker (i.e. seed point) for locating the pulmonary trunk is inside of the Pulmonary Artery then, it is considered as successfully locating the pulmonary trunk. One may then segment this structure until one reach the pulmonary valve, and designated all those tracked regions as a pulmonary artery (pulmonary trunk), to extract the whole structure of the pulmonary trunk. The main pulmonary arteries are proximal to the ascending aorta, and remain on its left. By using this knowledge, one may first create a rectangular search space adjacent to the left lateral side of ascending aorta, then apply the segmentation pipeline of the invention to the rectangular search space in order to detect and segment a region which resembles pulmonary artery morphology. Further, this region (possible pulmonary artery candidates) is tracked in the cranial to caudal direction to the conus arteriosus in order to reach level of the pulmonary valve (PV). (The conus arteriosus, also known as infundibulum, is a conical pouch formed from the upper and left angle of the right ventricle in the chordate heart, from which the pulmonary trunk arises, and it develops from the bulbus cordis.) The circularity of the segmented region is used to determine whether level of the PV has been reached. If the level of the PV was reached, then one may designate all those tracked regions as a pulmonary artery. Since pulmonary arteries obtained next to the ascending aorta between the level of carina trachea to the apical level of PV, these arteries are the main pulmonary arteries which is also named as the pulmonary trunk.

The segmentation pipeline consist of three steps: -Anisotropic diffusion filtering for image enhancement, -Eigenvalues of Hessian matrix and Canny edge detection filter for finding boundary of structures and

-Morphological operations, the gray scale segmentation and 2D region growing for grouping pixels in a one segment.

Anatomically, the pulmonary trunk appears next to the ascending aorta in the CT slice. Figure 11 shows a search space for the pulmonary trunk (PT), where the area delineated by the shaded rectangle in the upper right of the image is the look up area for the PT in the CT slice. Once one have located the ascending aorta, one can locate the pulmonary trunk by creating a rectangular search space adjacent to the ascending aorta based on the ascending aorta location in the I (: ,: ,z_aa+2) (Figure 11). The location and the area of the rectangular depends on the location and the size of the ascending aorta in the CT slice. The length of the rectangle was set to 66 pixel from empirical observations. However, the width (w, in pixel) of the rectangle is calculated by w = | eAot-eAob | +15 [Equation 15] where eAo_t is the topmost and eAob the bottom most pixel of the ascending aorta region. To ensure that the search space will cover the pulmonary trunk despite inter-individual anatomical differences, one can extend the width by 15 pixels. The pixel location (x and y coordinates) of the rectangle is calculated by: r(ymin -.ymax ^min ■. _max)= ((eAo_t - 15) eAob , (e4o/+5) : (6i4o/+70)) [Equation 16] where eAoi is the left most pixel of the ascending aorta region. One may then draw the rectangle 15 pixels from above from the top point of the ascending aorta and from 5 pixels left of the left most pixel. The use of 5 and 15 pixels were determined by empirical observation. The following methods may be used: Anisotropic diffusion filtering, Eigen values of Hessian matrix, Canny edge detection, and Morphological operations.

After the pulmonary trunk has been located/ detected, detection of the pulmonary valve is the next step. Detecting (S17) an apical level of pulmonary valve comprises: segmenting and tracking (S17a) the pulmonary trunk slice by slice from carina level to conus arteriosus in cranial to caudal direction, identifying (S17b) a region segmented in the last slices from step 17a, evaluating (S17c) sphericity of the segmented region using a first criterion, cirl, to identify the pulmonary valve, on condition that the pulmonary valve is not identified using the first criterion, evaluating (S17d) the sphericity of the segmented region using a second criterion, cir2, to identify the pulmonary valve, detecting (S17e) an apical level of the pulmonary valve as the last segmented component identified in the previous step S16b above.

It is difficult to determine the location of pulmonary valve in CTPA volume images. However, anatomically is known that pulmonary valve is inside of the circular object (in the axial plane) somewhere between the right ventricle and the pulmonary artery. Therefore, a circular/roundish proximal part of the pulmonary trunk can be regarded as the apical level of the pulmonary valve.

In the previous step the pulmonary trunk in the CT stack /(: ,: ,z_pt) has been located. One purpose here is to calculate the diameter of the pulmonary trunk in a given CT scan. Calculating the diameter of the pulmonary trunk on several slices gives more accurate and precise results instead of calculating in one slice. Thus, one may track the pulmonary artery starting from the pulmonary trunk at the level of the carina trachea up to the conus arteriosus (infundibulum) in order to reach the apical level of pulmonary valve in the cranial to caudal direction. Anatomically, if one track the upper areas of the pulmonary trunk in cranial to caudal direction, one can reach the conus arteriosus, where the apical level of pulmonary valve is located. Therefore, one first take a trackable region from the upper part of the pulmonary trunk (Figure 11), and then segment and track the pulmonary trunk slice by slice using the following methods: Gray scale segmentation, Anisotropic diffusion filtering, Eigen values of Hessian matrix, Canny edge detection, Morphological operations, and 2D region growing. Figure 12 shows a tracking area for the pulmonary trunk, where the area delineated by the grey dot in the circle is the sampling area for tracking the pulmonary trunk up to the level of the apical level of pulmonary valve in cranial to caudal direction.

In the previous step, the method segment and track the pulmonary artery starting from the pulmonary trunk at the level of the carina trachea up to the conus arteriosus (infundibulum) in order to reach the apical level of pulmonary valve. To know that the last region that was segmented is the apical level of pulmonary valve, one may check the circularity of the segmented region. To check circularity (or sphericity), consider the two criteria below.

The first criteria is: ciri = AnA/ P² [Equation 17] where A and P are the area and the perimeter of the component respectively. The output of the first criteria is between 0 and 1. A component whose ciri is 1 is actually a circular object, but since the pulmonary trunk has a complex morphology, ciri of the apical level of pulmonary valve will not be 1. It was empirically determined a ciri cut-off value of 0.75 to designate the component as an apical level of pulmonary valve. If the calculated ciri value is < 0.70, continue to segment and track the pulmonary artery until a circular object is found. However, if the calculated ciri value is between 0.70 and 0.75, the component may not be identified as an apical level of pulmonary valve. Accordingly, the apical level of pulmonary valve can be missed by the first criteria.

To find the apical level of pulmonary valve in case it was missed by the first criteria, a second criteria (cir2) is defined consisting of three sub steps. First, apply Hough transform to the segmented component in order to find vertical and horizontal lines inside of the object (Figure 13), where figure 13 shows the horizontal (image on the left) and the vertical (image on the right) lines that calculated by the Hough transform. Second, search for the longest vertical and horizontal lines by comparing line lengths. Finally, the absolute value difference of the longest vertical and the longest horizontal lines is calculated. If the calculated absolute value difference is < 20 pixels, then the component is designated as a circular object. Theoretically, the length of the major (horizontal) and the minor (vertical) axes of an ellipse is equal. Hence, empirically it was found that the cut-off value 20 is adequate to designate component as a circular object. So one can formulize circularity checking algorithm below:

After detecting the apical level of pulmonary valve, one may extract the apical level of pulmonary valve. The last segmented component by the previous step in the CT slice /(: ,z_pv) is the apical level of pulmonary valve.

After locating the apical level of pulmonary valve, in the next step one may determining (S18) a representative diameter of the pulmonary trunk by segmenting (S18a) the pulmonary trunk, and tracking (S18b) the segmented pulmonary trunk, slice by slice, from carina level to the apical level of pulmonary valve and back to carina level, wherein a diameter of the pulmonary trunk is measured in each slice; and calculating (S18c) a representative diameter of the pulmonary trunk as a mean of the measured diameters from each slice.

Optionally slices at the ends of the segmented pulmonary trunk may be removed when calculating the mean using heuristics. The Hough transform is applied to the all tracked pulmonary arteries (potentially pulmonary trunk since main pulmonary arteries are the pulmonary trunk) to calculate diameters. Roughly, the Hough Transform is used here for drawing a line (horizontal or vertical line with a given angle) between the opposite sides of a polygon. Here, the length of the line is the potential diameter. The diameter is measured in each slice on the upper part of the slice in axial view, the spine facing downwards, and the mean diameter calculated as the average + a standard deviation. Finally, the representative pulmonary trunk diameter is the mean of all diameters calculated by the Hough Transform.

Example Device and System Configurations

Turning now to Figure 14, which is a block diagram of a computing device of the present disclosure. The computing device is configured to implement all aspects of the methods described in relation to Figure 4, e.g. by using a software where the method of the invention is stored. The computing device 10 comprises a communication interface (i/f) 11 configured for communication e.g. with a server for obtaining CT data.

The computing device 10 comprises a controller, CTL, or a processing circuitry 12 that may be constituted by any suitable Central Processing Unit, CPU, microcontroller, Digital Signal Processor, DSP, etc. capable of executing computer program code. The computer program may be stored in a memory, MEM 13. The memory 13 can be any combination of a Read And write Memory, RAM, and a Read Only Memory, ROM. The memory 13 may also comprise persistent storage, which, for example, can be any single one or combination of magnetic memory, optical memory, or solid state memory or even remotely mounted memory.

The present disclosure provides a computer-aided detection, CADe, method performed in a computing device for locating and measuring large vessels and detecting ascending aorta hypertension and/or pulmonary trunk hypertension in a subject, the method comprising: performing the method of steps S10-S18 above, to determine a representative diameter of an ascending aorta and/or pulmonary trunk in a subject, comparing the determined representative diameter(s) with a preset threshold, determining that ascending aorta hypertension and/or pulmonary trunk hypertension is present in a subject if the determined respective representative diameter(s) is greater than the preset threshold. The pulmonary trunk and ascending aorta will have different thresholds for comparison in view of their respective diameters, and also the ratio between the diameter pulmonary (main pulmonary artery) and the ascending aorta may be used as a marker (threshold) to detect a condition in the subject. The thresholds set by the system may depend on the subject being scanned, where age, sex and size of the subject being scanned may be taken into account if applicable. Mostly, the thresholds are experimentally determined. In one example, a determined average diameter of the pulmonary arteries above 29 mm (threshold set to 29 mm) is indicative of pulmonary hypertension, which may be due to PE, heart disease, thrombo-emboli, etc. (applicable to adults of both genders), see e.g. Sanai, Shirin MD; Aronow, Wilbert S. MD; Ravipati, Gautham MD; Maguire, George P. MD; Belkin, Robert N. MD; Lehrman, Stuart G. MD Prediction of Moderate or Severe Pulmonary Hypertension by Main Pulmonary Artery Diameter and Main Pulmonary Artery Diameter/Ascending Aorta Diameter in Pulmonary Embolism, Cardiology in Review: September-October 2006 - Volume 14 - Issue 5 - p 213-214. In literature, there are many meta-analyses shows cut-off values (thresholds), see e.g. Pena, E., Dennie, C., Veinot, J., & Muniz, S. H. (2012), Pulmonary hypertension: how the radiologist can help, Radiographics, 32(1), 9-32, "The ratio of the main pulmonary arterial diameter to that of the ascending aorta is also greater than or equal to 1, another useful sign of pulmonary hypertension", or Edwards, P. D., Bull, R. K., & Coulden, R. (1998), CT measurement of main pulmonary artery diameter, The British journal of radiology, 71(850), 1018-1020. Thus, the skilled person would know which threshold to set for a respective person when carrying out the method on a computing device. A summary of potential marker values (thresholds) that may be used is summarized in table 2 below.

Table 2

In table 2, the mPAd is the diameter of the main pulmonary artery (pulmonary trunk) and AAd is the diameter of the ascending aorta, COPD is chronic obstructive pulmonary disease, HFpEF is heart failure with preserved ejection fraction, and RVD is right ventricular dysfunction.

A presence of ascending aorta hypertension and/or pulmonary trunk hypertension in the subject indicates presence of a medical disorder is the subject, the medical disorder being one or more of pulmonary embolism, heart disease, thrombo-embolism, cryptogenic organizing pneumonia, scleroderma, chronic obstructive pulmonary disease, heart failure with preserved ejection fraction and right ventricular dysfunction.

Thus, it is provided a computing device (10) comprising a memory (13) for storing instructions, and processing circuitry (12) for executing the instructions, wherein the processing circuitry (12) is configured to perform the methods above and below.

The present disclosure further provides a computer-aided detection system (200) comprising a computing device (10), configured for performing a computer-aided detection, CADe, method of locating and measuring large vessels in a subject undergoing a computed tomography, CT, scan, the computing device (10) comprising: a memory (13) memory for storing instructions and processing circuitry (12) configured to cause the computing device (10): to obtain CT data (20) from a consecutive CT scan of thorax of the subject, the CT data (20) comprising an image stack forming a CT volume of images, to prepare the CT data (20) for segmentation tasks, to locate carina trachea in the CT volume, to detect descending aorta at carina level, to detect ascending aorta, to detect pulmonary trunk, and to detect an apical level of pulmonary valve.

The system may be configured to perform any one of the methods above and below.

It further provided a computer-aided detection system (200) configured for determining a representative diameter of a pulmonary trunk in a subject, the system comprising: a processor (14) configured to determine a representative diameter of a pulmonary trunk in a CT volume of images from a CT scan of the subject by tracking a segmented pulmonary trunk from carina level to pulmonary valve and back to carina level, wherein a diameter of the pulmonary trunk is measured in each slice, and a representative diameter of the pulmonary trunk is calculated as a mean of the measured diameters from each slice.

According to some aspects, the disclosure relates to a computer program comprising computer program code which, when executed, causes a computing device to execute the methods described above and below. According to some aspects the disclosure pertains to a computer program product or a computer readable medium holding said computer program. The processing circuitry may further comprise both a memory 13 storing a computer program and a processor 14, the processor being configured to carry out the method of the computer program. According to some aspects is provided a carrier containing any one of the computer programs mentioned above, wherein the carrier is one of an electronic signal, optical signal, radio signal, or computer readable storage medium.

The content of this disclosure thus enables automatic location and measuring of large vessels in a CT scan volume of images of a subject, which may be used for automatic detection of a medical condition in the subject.

Aspects of the disclosure are described with reference to the drawings, e.g., block diagrams and/or flowcharts. It is understood that several entities in the drawings, e.g., blocks of the block diagrams, and also combinations of entities in the drawings, can be implemented by computer program instructions, which instructions can be stored in a computer-readable memory, and also loaded onto a computer or other programmable data processing apparatus. Such computer program instructions can be provided to a processor of a general purpose computer, a special purpose computer and/or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer and/or other programmable data processing apparatus, create means for implementing the functions/acts specified in the block diagrams and/or flowchart block or blocks.

In the drawings and specification, there have been disclosed exemplary aspects of the disclosure. However, many variations and modifications can be made to these aspects without substantially departing from the principles of the present disclosure. Thus, the disclosure should be regarded as illustrative rather than restrictive, and not as being limited to the particular aspects discussed above. Accordingly, although specific terms are employed, they are used in a generic and descriptive sense only and not for purposes of limitation.

The description of the example embodiments provided herein have been presented for purposes of illustration. The description is not intended to be exhaustive or to limit example embodiments to the precise form disclosed, and modifications and variations are possible in light of the above teachings or may be acquired from practice of various alternatives to the provided embodiments. The examples discussed herein were chosen and described in order to explain the principles and the nature of various example embodiments and its practical application to enable one skilled in the art to utilize the example embodiments in various manners and with various modifications as are suited to the particular use contemplated. The features of the embodiments described herein may be combined in all possible combinations of methods, computing devices, modules, systems, and computer program products. It should be appreciated that the example embodiments presented herein may be practiced in any combination with each other.

It should be noted that the word "comprising" does not necessarily exclude the presence of other elements or steps than those listed and the words "a" or "an" preceding an element do not exclude the presence of a plurality of such elements. It should further be noted that any reference signs do not limit the scope of the claims, that the example embodiments may be implemented at least in part by means of both hardware and software, and that several "devices" may be represented by the same item of hardware.

The various example embodiments described herein are described in the general context of method steps or processes, which may be implemented in one aspect by a computer program product, embodied in a computer-readable medium, including computerexecutable instructions, such as program code, executed by computers in networked environments. A computer-readable medium may include removable and non-removable storage devices including, but not limited to, Read Only Memory (ROM), Random Access Memory (RAM), compact discs (CDs), digital versatile discs (DVD), etc. Generally, program modules may include routines, programs, objects, components, data structures, etc. that performs particular tasks or implement particular abstract data types. Computer-executable instructions, associated data structures, and program modules represent examples of program code for executing steps of the methods disclosed herein. The particular sequence of such executable instructions or associated data structures represents examples of corresponding acts for implementing the functions described in such steps or processes.

EXAMPLES

EXAMPLE 1, AUTOMATIC DETECTION OF LARGE VESSELS AND MEASUREMENTS THEREOF IN COMPARISON TO MANUAL LABOR

A single-center retrospective study included 700 non-ECG-gated CTPA examinations from 652 patients (median age, 72 years and interquartile range 18 years; age range 16-100 years; 353 women) performed between 2014 and 2018. A training set of 180 examinations was used to develop a fully automated deterministic algorithm, and the final model was tested using 520 examinations. The detection and segmentation performance was qualitatively validated by visual evaluation.

Dataset

The dataset consists of 700 retrospective non-ECG-gated CT pulmonary angiography examinations performed at a single institution (Nykbping Hospital, Nykbping, Sweden) between 2014 and 2018. 383 CTPA examinations from 353 women (age range 16-97 years; median age 73 years; interquartile range 20 years) and 317 from 299 men (age range 19-100 years; median age 71 years; interquartile range 15 years). The CTPAs were randomly assigned to training (n=180) and test sets (n=520). The training set was used for algorithm development while the test set was prepared for unbiased evaluation of the final algorithm.

CT Image Acquisition

Non-ECG-gated CTPA was performed with 5 different multidetector-row CT scanners (Brilliance 64, Ingenuity Core and Ingenuity CT, Philips Medical Systems, Eindhoven, the Netherlands; Lightspeed VCT, General Electric (GE) Healthcare Systems, Waukesha, Wl, USA; Somatom Definition Flash, Siemens Healthcare, Erlangen, Germany). Examinations were performed after injection of intravenous (IV) contrast (Omnipaque 350 mgl/ml, GE Healthcare Systems, Waukesha, Wl, USA) and saline. The CT image acquisition technique varied by manufacturer with most frequent slice thickness of 0.625 mm (0.625 mm - 2.0 mm), pixel spacing of 0.7 mm (0.59 mm - 0.98 mm), and voltage of 100 kV (80 kV - 120 kV). A secondary axial reformat with 2.0 mm slice thickness was performed on all CTPA examinations.

Manual Measurements

The CTPA examinations were exported from the Picture Archiving and Communication System (PACS, Sectra AB, Linkbping, Sweden) system in DICOM format. The examinations were reviewed and annotated using the RadiAnt DICOM Viewer software (Medixant) by a senior radiologist (TF) with 15 years of experience. 150 of the CTPAs were first reviewed and annotated by a radiology resident (DT) with 5 years of experience in general diagnostic radiology and then double read by TF. For all CTPAs the final annotation was decided by TF. Scrolling in the 2 mm axial image stack, the image which optimally presented the pulmonary trunk (PT) was identified and in this image the following manual measurements were performed (Figure 15): The diameter of the PT and ascending aorta (AAo), IV contrast concentration in PT (mean value of HU in 2 cm² circular region of interest, ROI), and image noise (SD of HU in a 1 cm² circular ROI in the descending aorta, DAo). Figure 15 show a demonstration of manual measurements by a radiologist. Figure 15A. IV contrast concentration in PT is calculated by taking the mean value of HU in 2 cm² circular region of interest (arrow). Image noise is calculated by taking the SD of HU in a 1 cm² circular region of interest in the descending aorta (arrow head). Figure 15B. The diameter of the PT (arrow) and ascending aorta (arrow head) are calculated by scrolling in the 2 mm axial image stack, the image which optimally presented the pulmonary trunk was identified. These measurements were used as ground truth (Radiological characteristics of 700 CTPA examinations used in CADe system training and testing). For each CTPA examination the radiologist also scored five different image quality parameters affecting the evaluation for PE, namely motion artifacts, streak artifacts, IV contrast concentration in pulmonary trunk, parenchymal disease and image noise. A total score was calculated, and the result determined the overall CTPA examination quality as good, acceptable or inferior (Figure 16). Figure 16 shows a demonstration of CT pulmonary angiography (CTPA) image quality examples at three different levels with the same windowing parameters (width = 600 HU, level = 100 HU).

CADe system

To automatically measure morphometric and geometric parameters of mediastinum structures (AAo, DAo, PT), we combined image processing and image analysis techniques to detect and segment compartments. The developed CADe system consists of two main steps, pre-processing and segmentation chain (Figure 17). Figure 17 illustrates a flowchart of the CADe algorithm with the parts A, Bl and B2, C1-C6, D1-D6, E1-E4, and F1-F8, as mentioned below.

Pre-Processing

After acquiring the 2 mm axial CT images, every voxel was converted to Hounsfield Units (HU) and the direction of scanning was determined based on information in the DICOM header. As patient orientation may vary along the cranial to caudal direction because of scoliosis, movement or position during the examination, or other reasons, the CT exam was aligned with respect to the X axis of the axial plane in the cranial as well as caudal part of the examination (Figure 17, Bl and B2). Segmentation Pipeline

First, the system automatically located the three-dimensional positions of two anatomical landmarks, the carina of trachea and an apical level of the pulmonary valve (PV), to find seed points of vascular structures in the mediastinum (AAo and DAo, PT). Second, seed points for segmentation of mediastinal structures were automatically placed using a heuristic approach. Third, the structures were segmented by applying the following image processing techniques; image enhancement, edge detection, gray scale segmentation, and 2D region growing methods. Finally, measurements were done on the segmented regions.

Locating the trachea and airways

The trachea can be automatically detected based on two simple anatomical features. The first of the air-filled structures encountered in the thoracic cavity when moving cranially to caudally is probably the trachea or one of the three air-filled structures in the thoracic cavity, of which the middle is likely the trachea and the other two correspond to the right and left lung. However, in the planes of some CTPA examinations, there are craniofacial structures or upper airway structures that resemble the trachea, artefacts such as beam hardening, or adjoining of the trachea to the right lung cranial to the carina which complicates accurate identification of the trachea. To overcome this problem and to correctly locate the trachea, the system extracts two adjoined volumes of interest from the superior part of the thoracic cavity, wherein trachea candidates are generated and assessed separately (Figure 17, C4 and C5). In these two volumes, the air-filled structures are identified by thresholding, flood-fill operations and connected component analyses. Trachea candidates are then joined across the two volumes and connected component analyses are applied. The longest of the candidates also having diameter and volume within an empirically determined range is selected as the trachea. As the trachea has been detected, it can be tracked cranially to caudally slice by slice to the bifurcation point where the left and right main bronchi can be found as two distinct segments (Figure 17, C6). We then designate the CT slice, where the distance between the left and right main bronchus is >0.75 cm, as the level of the carina of the trachea. Thus, the trachea and the carina of the trachea can be located automatically in 3D. Locating major vascular structures of the mediastinum.

We observed that the DAo can be easily detected around the level of carina trachea since the DAo is always located posterior to the left main bronchus anatomically and its morphological appearance (IV contrast concentration and circularity) is homogeneous around the level of carina trachea because of CTPA contrast injection protocols and human anatomy. Therefore, we generate eight artificial rays that search spaces posterior to the left main bronchus at the level of carina trachea (Figure 17, DI). However, a second set of eight rays is needed as the DAo is occasionally shifted toward the spinal canal. The mass center of the DAo can then be located to the right of the mass center of the left main bronchus, which causes the first set of rays to miss the DAo. Therefore, two sets of ray search spaces are sequentially used and one or more rays passing a sufficient number of connected pixels to correspond to the contrast filled aorta indicate that the DAo has been located (Figure 17, D3). Once the DAo has been detected, the segmentation process was done as described above. The area (1 cm²) enclosed by a circle at the mass center of DAo was used to calculate the average HU density of the CTPA examination as a metric of contrast filling (Figure 17, D6).

Once the DAo has been segmented, the AAo can be easily detected by tracking and comparing the segmented DAo regions slice by slice in the caudal to cranial direction to find geometrical features characteristic of the AAo. Morphological changes indicate the level of the aortic arch (Figure 17, El). By tracking the upper part of the aortic arch in the cranial to caudal direction slice by slice the first identified circular object was designated as AAo (Figure 17, E2). Next, we segment the AAo between the levels of the aortic arch and the carina of trachea (Figure 17, E3) and calculate the diameter of AAo as the mean of AAo diameters in these planes.

The pulmonary trunk (PT) is adjacent to the ascending aorta and remain on its left. We therefore first created a rectangular search space adjacent to the left lateral side of AAo (Figure 17, F2) and applied our segmentation pipeline to said search space to segment the pulmonary trunk. Next, the pulmonary trunk was tracked in the cranial-to-caudal direction in order to reach the apical level of the pulmonary valve (PV). The circularity of the segmented region was used to determine whether the apical level of the PV had been reached (Figure 17, F4). The tracked pulmonary trunk diameter was calculated by Hough transform (Figure 17, F4) as the mean of PT diameters in these planes. Taken together, the DAo, AAo and PT were automatically detected in 3D and their average diameters and contrast levels obtained.

Statistical Analysis

Statistical analysis was performed using MATLAB R2019b and Microsoft Office Excel. Bland- Altman and scatter analyses were used to compare agreement and relationship, respectively, between automated and manual measurements. A p-value of <.05 was defined as statistically significant and Spearman's correlation coefficient was used to evaluate the agreement between automated measurement and manual measurement.

Results and Conclusion

Measuring morphometric parameters, such as Hounsfield unit radiodensity, and geometric parameters such as diameters of mediastinal vessels in CT imagery is time consuming for the radiologist but can aid diagnosis of cardiovascular conditions like pulmonary embolism and pulmonary hypertension. The developed algorithm accurately and automatically measured the diameter of the pulmonary trunk and ascending aorta, the IV contrast concentration in the pulmonary trunk and image noise in less than 10 s per CTPA examination.

The algorithm correctly located the carina of trachea in 89% of examinations, while the apical level of the pulmonary valve was detected in 83%. The descending aorta (DAo), the ascending aorta (AAo) and the pulmonary trunk (PT) was detected in 89%, 87%, and 84% of examinations, respectively. For correctly detected vessels, qualitative evaluation revealed successful segmentation in 100% for the DAo, 99% for the AAo, and 98% for the PT. The automatic measurements correlated well with those of the radiologist (Spearman's r = 0.89 for image noise, 0.99 for IV contrast in the PT, 0.92 for AAo diameter, and 0.71 for PT diameter). The CTPA image quality slightly affected detection of landmarks and mediastinal vessel structures, but did not affect the computational performance of the algorithm. Thus, the fully automated algorithm accurately detected, segmented, and measured mediastinal vessel structures in large number of unselected CTPA examinations with an adequate representation of common artifacts and medical conditions. EXAMPLE 2, EVALUATION OF SEGMENTATION RESULTS

The main goal of the study on Example 1 is to compare radiologist measurements versus developed algorithm measurements. We use numerical data (ground truth) to make a quantitative comparison of radiologist measurements versus developed algorithm measurements.

To make a valid measurement, good image segmentation is needed. The good measurement results are enough to estimate that the segmentation results are also good. However, to test the strength of the developed algorithm and to increase the quality of the study, evaluation of the image segmentation is required. The common way to evaluate image segmentation is using quantitative evaluation methods such as Dice Coefficient, Jaccard's Index, or Pixel accuracy. But all these methods require pixel-wise annotations, and it is time-consuming. Making a pixel-wise annotation of more than one class in a large data set requires overwhelming work. Since we have a large dataset (n=700) and multi-object (n=4) segmentation task, a remaining option for us is to use the qualitative evaluation method. Here, we defined the evaluation criteria to evaluate image segmentation results.

We presented segmentation results in two ways, using a mask and a boundary box. The mask or boundary box with a tag was inserted on a CT image by our developed algorithm that represents which organ was segmented. The mask contains a set of pixels and the boundary box is a rectangle that covers an object.

For Evaluating segmentation results checking the mask or the boundary box is adequate. Here, we encouraged the radiologists to evaluate with the boundary box for fast and easy evaluation. The radiologist's mission is here to check if the boundary box covers the desired organ or not. However, if the radiologists cannot decide that the segmentation is passed or failed by evaluating with the boundary box (for instance the boundary box covers also other organs/tissues), the radiologists can go deeper by evaluating with the mask.

In this study, we have 4 segmentation tasks, namely:

Task 1: Segmentation of the Trachea,

Task 2: Segmentation of the Descending Aorta, Task 3: Segmentation of the Ascending Aorta,

Task 4: Segmentation of the Pulmonary Trunk.

Each task has its own evaluation criteria depending on aim of the segmentation. The radiologists can find in the following sections.

Task 1: Segmentation of the Trachea

Evaluation Criteria: We have a 3D-Segmentation of the trachea from superior slices of the CT examinations to the bifurcation point. But we only evaluate one of the segmented trachea areas in 3D-segmentation volume.

• If the boundary box covers the majority (estimated 50% or more of the whole area) of trachea then, it can be labeled as passed otherwise it must be labeled as failed.

• If the boundary box over segmented (the box is clearly larger than the trachea) the trachea then, it must be labeled as failed.

Example: The boundary box covers the whole trachea, in that case, the segmentation task can be labeled as passed (Figure 18). Figure 18 shows an example of successful segmentation. In the upper image a), a boundary box is used, and it can be clearly seen that the boundary box covers the whole trachea. In the image b), the mask is used, and it can be clearly seen that the mask overlaps with only and whole trachea.

Task 2: Segmentation of the Descending Aorta

Evaluation Criteria:

• If the boundary box covers the majority (estimated 75% or more of the whole area) of descending aorta then, it can be labeled as passed otherwise it must be labeled as failed.

• If the boundary box over segmented (the box is clearly larger than the descending aorta) the descending aorta then, it must be labeled as failed.

Example: The boundary box covers the whole descending aorta, in that case, the segmentation task can be labeled as passed (Figure 19). Figure 19 is an example of successful segmentation. In the image a), a boundary box is used, and it can be clearly seen that the boundary box covers the whole descending aorta. In the image b), the mask is used, and it can be seen that the mask overlaps with the whole and only the descending aorta (mask is circular dark grey dot within the black circle).

Task 3: Segmentation of the Ascending Aorta

Evaluation Criteria:

We have a 3D-Segmentation of the ascending aorta. But we only evaluate one of the segmented areas in 3D-segmentation volume.

• If the boundary box covers the majority (estimated 75% or more of the whole area) of the ascending aorta then, it can be labeled as passed, otherwise it must be labeled as failed.

• If the boundary box over segmented (the box is clearly larger than the ascending aorta) the ascending aorta then, it must be labeled as failed.

Example: The boundary box covers the whole ascending aorta, in that case, the segmentation task can be labeled as passed (Figure 20). In the image, a boundary box is used, and it can be clearly seen that the boundary box covers the whole ascending aorta.

Task 4: Segmentation of the Pulmonary Trunk

Evaluation Criteria:

We have a 3D-Segmentation of the main pulmonary artery. But we only evaluate one of the segmented areas in 3D-segmentation volume.

• If the boundary box clearly covers the upper (y-axis, in the axial-plane) part of the main pulmonary trunk, then it can be labeled as passed, otherwise it must be labeled as failed.

• If the radiologists cannot decide that the segmentation is passed or failed by evaluating with the boundary box (for instance the boundary box covers also other organs/tissues), the radiologists can check the mask by following criteria: o If the mask covers the majority (estimated 75% or more of the whole area) of the pulmonary artery then, it can be labeled as passed, otherwise it must be labeled as failed. o If the mask over segmented (the mask is clearly larger than the pulmonary artery or includes other organs/tissues) the pulmonary artery, then it must be labeled as failed.

Example: The boundary box covers the whole main pulmonary artery (pulmonary trunk), in that case, the segmentation task can be labeled as passed (Figure 21). Figure 21 is an example of successful segmentation. In the image, a boundary box is used, and it can be clearly seen that the boundary box covers the whole main pulmonary artery.

EXAMPLE 3, EVALUATION OF THE DETECTION RESULTS

The marker which was inserted on a CT image by our developed algorithm represents which organ or anatomical landmark we detected. This marker consists of a single pixel also known as seed point. We have resized the marker (becoming filled circle) so that this point can be seen more easily on the image. If the marker hits one of the voxels of the desired organ, then the detection task is accomplished. Here, the radiologist's mission is to check if the marker is in the desired organ or not. If it is in the desired organ, then the detection task is labeled as passed, otherwise failed.

In this study, we have 7 detection tasks, namely:

Task 1: Detection of the Trachea,

Task 2: Detection of the Bifurcation Point of Airways

Task 3: Detection of the Carina Level,

Task 4: Detection of the Descending Aorta,

Task 5: Detection of the Ascending Aorta,

Task 6: Detection of the Apical Level of Pulmonary Valve,

Task 7: Detection of the Pulmonary Trunk.

Each task has its own evaluation criteria depending on aim of the detection, a specified in the following sections, but there is a general criterion for failed cases, which is described below. Task 1: Detection of the Trachea

Evaluation Criteria: If the marker is inside of the Trachea then, it can be labeled as passed otherwise it must be labeled as failed.

Task 2: Detection of the Bifurcation Point of Airways

Evaluation Criteria: If the two markers are inside of the left and the right main bronchus then, it can be labeled as passed otherwise it must be labeled as failed.

Task 3: Detection of the Carina Level

Evaluation Criteria: According our observation, if the distance between the left main bronchus and the right main bronchus is bigger than 0.7 cm and smaller than 1.7 cm then we can easily be able to detect the descending aorta. So, in the image you can be able to see the distance, If the distance is between 0.7 cm and 1.7 cm then, it can be labeled as passed otherwise it must be labeled as failed.

Task 4: Detection of the Descending Aorta

Evaluation Criteria: If the marker is inside of the Descending Aorta then, it can be labeled as passed otherwise it must be labeled as failed.

Task 5: Detection of the Ascending Aorta

Evaluation Criteria: If the marker is inside of the Ascending Aorta then, it can be labeled as passed otherwise it must be labeled as failed.

Task 6: Detection of the Apical Level of Pulmonary Valve

What is it that we call the Apical Level of the Pulmonary Valve?

It is hard to determine the location of pulmonary valve in CTPA volume images. However, anatomically we know that pulmonary valve is inside of the circular object (in the axial plane) somewhere between the right ventricle and the pulmonary artery.

Assume that, we find an object that meets above criteria. Therefore, a circular object that does not touch (separated) with the left or the right pulmonary arteries but belongs to the pulmonary artery system can be named as Apical Level of Pulmonary Valve.

Evaluation Criteria: If the marker is inside a circular object which belongs to pulmonary artery system then, it can be labeled as passed otherwise it must be labeled as failed.

Task 7: Detection of the Pulmonary Trunk

Evaluation Criteria: If the marker is inside of the Pulmonary Artery then, it can be labeled as passed, otherwise it must be labeled as failed.

EMBODIMENTS

1. A computer-aided detection, CADe, method performed in a computing device for locating and measuring large vessels in a subject undergoing a computed tomography, CT, scan, the method comprising: obtaining (S10) CT data from a consecutive CT scan of thorax of the subject, the CT data comprising an image stack of slices forming a CT volume of images; preparing (Sil) the CT data for segmentation tasks, to ensure that the orientation of the subject in view of the CT scan is within an acceptable range; locating (S12) carina trachea in the CT volume; detecting (S13) descending aorta at carina level; detecting (S14) ascending aorta; detecting (S16) pulmonary trunk; and detecting (S17) an apical level of pulmonary valve.

2. The method of embodiment 1, further comprising: determining (S15) a representative diameter of the ascending aorta.

3. The method of embodiment 1, further comprising: determining (S18) a representative diameter of the pulmonary trunk. The method of embodiment 2, wherein determining (S15) a representative diameter of the ascending aorta comprises: segmenting (S15a) the ascending aorta; tracking (S15b) the ascending aorta slice by slice from the first located slice to the carina trachea, wherein a diameter of the segmented ascending aorta is measured in every slice; and calculating (S15c) a representative diameter of the ascending aorta as a mean of the measured diameters from each slice. The method of embodiment 3, wherein determining (S18) a representative diameter of a pulmonary trunk comprises: segmenting (S18a) the pulmonary trunk, and tracking (S18b) the segmented pulmonary trunk, slice by slice, from carina level to apical level of pulmonary valve and back to carina level, wherein a diameter of the pulmonary trunk is measured in each slice; and calculating (S18c) a representative diameter of the pulmonary trunk as a mean of the measured diameters from each slice. The method of any one of embodiments 1-5, wherein obtaining (S10) the CT data comprises: retrieving, using the computing device, stored image matrices as CT data from a server to obtain the CT data, wherein the CT data is produced by: performing a CT scan of the subject using a CT scanner to retrieve a number of slices; sending the retrieved slices to a computing device for reconstruction; turning each slice into an image matrix using a software in the computing device; and reformatting and storing the image matrices under DICOM standards as CT data in one or more servers. The method of any one of embodiments 1-3, wherein preparing (Sil) the dataset for segmentation tasks comprises: calculating (Sila) linear scale value; obtaining (Sllb) information indicating a scanning direction of the CT scan; determining (sllc) an orientation of the subject, wherein determining the orientation of the subject comprises calculating the subject's orientation in an image of the image stack as an angle a between a major axis of the image and an x-axis of an image plane; and on condition that a is larger than 17 °, performing (Slid) an additional calculation of the subject's orientation based on caudal slices of the CT scan to attain an updated angle a; and correcting (Slle) scan curvature by rotating the CT scan according to a.

8. The method of any one of embodiments 1-7, wherein locating (S12) the carina trachea in the CT volume by locating trachea at carina level comprises: locating (S12a) trachea in the CT volume; detecting (S12b) any tracheal intubation present; segmenting (S12c) airways of trachea region; and identifying (S12d) the optimal carina location for the trachea at carina level as the carina trachea.

9. The method of embodiment 8, wherein locating (S12a) the trachea in the CT volume comprises: downscaling the CT volume to half its size; determining a first volume of interest, VOI1; segmenting air areas in VOI1; generating a pool of potential trachea candidates, TCP1, in VOI1; and on condition that there is one 3-dimensional, 3D, component in the TCP1, locating the trachea as the 3D component; on condition that there is more than one 3D component in the TCP1, determining a second volume of interest, VOI2; segmenting air areas in VOI2; generating a second pool of potential trachea candidates, TCP2, in VOI2; combining TCP1 and TCP2 to find common 3D components; locating the trachea as the 3D component spanning the largest total number of CT images of the combined TCP1 and TCP2; and selecting a slice from the CT volume of images where the trachea is located as a seed point S_t for a subsequent extraction of the trachea, wherein detecting (S12b) any tracheal intubation present comprises: calculating a center point of the trachea in a binary image from the image stack; performing image thresholding on the binary image; excluding areas in the binary image having less than 15 and above 100 pixels, i.e. areas being smaller than 7,56 mm² and larger than 50,4 mm²; applying 2D connected component analysis to the binary image; calculating center points of all distinct 2D connected components in the binary image; calculating a Euclidean distance of component center points to the center point of the trachea to identify the closest components to the trachea center point; identifying any center point of the 2D connected components having a distance of less than 11 pixels, corresponding to an area of 5,55 mm² or less, as a potential endotracheal tube; and on condition that a tracheal intubation is present, performing 3D region growing to the CT data starting at the seed point S_t to extract the tracheal intubation; and identifying a slice in the image stack of the extracted tracheal intubation in cranial to caudal direction as the last image where the tracheal intubation exists, and selecting said slice as seed point S_t for a subsequent extraction of the trachea, wherein segmenting (S12c) the airways comprises: thresholding the selected slice to obtain a binary image; and segmenting and tracking the trachea region up to bifurcation point in cranial to caudal direction using the 2-dimensional, 2D, growing method using the binary image and the seed point S_t as inputs, wherein the selected slice and seed point S_t is taken from step 12a if no tracheal intubation exists, and from 12b if tracheal intubation exists, wherein identifying (S12d) the optimal carina location as the carina trachea comprises: identifying the optimal carina location by tracking the trachea slice by slice from the bifurcation point in cranial to caudal direction, calculating the distance between left main bronchus and right main bronchus for each respective slice until reaching a slice, Soc, where the calculated distance is equal to or greater than 0,75 cm, and selecting said slice Socas the optimal carina location; and identifying the selected optimal carina location as the carina trachea. The method of any one of embodiments 1-9, wherein detecting (S13) descending aorta at carina level comprises: creating a first search space originating from an artificial line that cuts the x- axis of the slice at an angle of 22,5 -75 ° downwards clockwise from the center of left main bronchus comprising lines having an angle of angle larger than 2,5° and smaller than 15° between each other forming rays from center of origin; creating a second search space originating from an artificial line that cuts the x-axis of the slice at an angle of 60 -120 ° downwards clockwise from the center of the left main bronchus comprising lines having an angle larger than 2,5° and smaller than 15° between each other forming rays from center of origin; identifying if any one of the rays pass through the descending aorta by going ray by ray and: collecting all pixels on the ray; thresholding said pixels over - 1 HU to acquire binary image; applying 2D connected component analysis to identify largest connected pixels on the ray; calculating the total number of pixels on the segmented component; and if the total number of pixels is greater than 15, i.e. the area is larger than 7,5 mm²; calculating mean x and y coordinates of those pixels, separately, in order to reduce to one data point, and identifying said data point as the descending aorta at carina level, PdesAo. The method of embodiment 10, wherein the lines of the first and second search spaces have an angle larger than 5° and smaller than 10° between each other, or an angle larger than 7° and smaller than 8° between each other, or an angle of 7,5° between each other. The method of any one of embodiments 1-11, wherein detecting (S14) ascending aorta comprises: detecting (S14a) aortic arch by tracking the descending aorta up to the aortic arch in caudal to cranial direction; detecting (S14b) the ascending aorta by tracking the upper part of the aortic arch to the ascending aorta in cranial to caudal direction to locate a first slice of the ascending aorta. The method of any one of embodiments 1-11, wherein detecting (S16) pulmonary trunk comprises: creating (S16a) a rectangular search space adjacent to the left lateral side of the ascending aorta location in the CT image stack covering the candidate of the pulmonary trunk, the pixel location of the rectangular being calculated as (eAot - 15) eAob, (E4OZ +5) : eAoi +70)), where eAo_t is the top most and eAob the bottom most pixel of the ascending aorta region; and detecting (S16b) the pulmonary trunk within said rectangular search space. The method of any one of embodiments 1-13, wherein detecting (S17) an apical level of pulmonary valve comprises: segmenting and tracking (S17a) the pulmonary trunk slice by slice from carina level to conus arteriosus in cranial to caudal direction; identifying (S17b) a region segmented in the last slices from step 17a; evaluating (S17c) sphericity of the segmented region using a first criterion, ciri, to identify the apical level of pulmonary valve; on condition that the apical level of pulmonary valve is not identified using the first criterion, evaluating (S17d) the sphericity of the segmented region using a second criterion, cir?, to identify the apical level of pulmonary valve, detecting (S17e) the apical level of the pulmonary valve as the last segmented component identified in step S16b.

15. The method of any one of the previous embodiments, wherein the image stack of consecutive CT scans are computed tomography pulmonary angiography (CTPA) volume images.

16. A computer-aided detection, CADe, method performed in a computing device for locating and measuring large vessels and detecting ascending aorta hypertension and/or pulmonary trunk hypertension in a subject, the method comprising: performing the method of any one of embodiments 2-15 to determine a representative diameter of an ascending aorta and/or pulmonary trunk in a subject; comparing the determined representative diameter(s) with a preset threshold; determining that ascending aorta hypertension and/or pulmonary trunk hypertension is present in a subject if the determined respective representative diameter(s) is greater than the preset threshold.

17. The method of embodiment 16, wherein a presence of ascending aorta hypertension and/or pulmonary trunk hypertension in the subject indicates presence of a medical disorder in the subject, the medical disorder being one or more of pulmonary embolism, heart disease, thrombo-embolism, cryptogenic organizing pneumonia and scleroderma.

18. A computing device (10) comprising a memory (13) for storing instructions, and processing circuitry (12) for executing the instructions, wherein the processing circuitry (12) is configured to perform the method of any one of embodiments 1 to 15.

19. A computer-aided detection system (200) comprising a computing device (10), configured for performing a computer-aided detection, CADe, method of locating and measuring large vessels in a subject undergoing a computed tomography, CT, scan, the computing device (10) comprising: a memory (13) for storing instructions and processing circuitry (12) configured to cause the computing device (10): to obtain CT data (20) from a consecutive CT scan of thorax of the subject, the CT data (20) comprising an image stack forming a CT volume of images; to prepare the CT data (20) for segmentation tasks; to locate carina trachea in the CT volume; to detect descending aorta at carina level; to detect ascending aorta; to detect pulmonary trunk; and to detect an apical level of pulmonary valve.

20. The computer-aided detection system of embodiment 19, configured to perform the method of any of embodiments 1-15.

21. A computer-aided detection system (200) configured for determining a representative diameter of a pulmonary trunk in a subject, the system comprising: a processor (14) configured to determine a representative diameter of a pulmonary trunk in a CT volume of images from a CT scan of the subject by tracking a segmented pulmonary trunk from carina level to pulmonary valve and back to carina level, wherein a diameter of the pulmonary trunk is measured in each slice, and a representative diameter of the pulmonary trunk is calculated as a mean of the measured diameters from each slice.

22. A computer program comprising computer program code which, when executed in a computing device, causes the computing device to execute the methods according to any of the embodiments 1-15.

23. A carrier containing the computer program of embodiment 22, wherein the carrier is one of an electronic signal, optical signal, radio signal, or computer readable storage medium. 24. The method, computing device, or computer-aided detection system according to any one of embodiments 1-21 above, wherein the CT is a contrast enhanced CT.

5

Claims

1. A computer-aided detection, CADe, method performed in a computing device for locating and measuring large vessels in a subject undergoing a computed tomography, CT, scan , the method comprising: obtaining (S10) CT data from a consecutive CT scan of thorax of the subject, the CT data comprising an image stack of slices forming a CT volume of images; preparing (Sil) the CT data for segmentation tasks; locating (S12) carina trachea in the CT volume; detecting (S13) descending aorta at carina level; detecting (S14) ascending aorta; detecting (S16) pulmonary trunk; and detecting (S17) an apical level of pulmonary valve.

2. The method of claim 1 further comprising: determining (S15) a representative diameter of the ascending aorta.

3. The method of claim 1 further comprising: determining (S18) a representative diameter of the pulmonary trunk.

4. The method of claim 2, wherein determining (S15) a representative diameter of the ascending aorta comprises: segmenting (S15a) the ascending aorta; tracking (S15b) the ascending aorta slice by slice from the first located slice to the carina trachea, wherein a diameter of the segmented ascending aorta is measured in every slice; and calculating (S15c) a representative diameter of the ascending aorta as a mean of the measured diameters from each slice.

5. The method of claim 3, wherein determining (S18) a representative diameter of a pulmonary trunk comprises: segmenting (S18a) the pulmonary trunk, and tracking (S18b) the segmented pulmonary trunk, slice by slice, from carina level to an apical level of pulmonary valve and back to carina level, wherein a diameter of the pulmonary trunk is measured in each slice; and calculating (S18c) a representative diameter of the pulmonary trunk as a mean of the measured diameters from each slice.

6. The method of any one of claims 1-5, wherein locating (S12) the carina trachea in the CT volume by locating trachea at carina level comprises: locating (S12a) trachea in the CT volume; detecting (S12b) any tracheal intubation present; segmenting (S12c) airways of trachea region; and identifying (S12d) the optimal carina location for the trachea at carina level as the carina trachea.

7. The method of claim 6, wherein locating (S12a) the trachea in the CT volume comprises: downscaling the CT volume to half its size; determining a first volume of interest, VOI1; segmenting air areas in VOI1; generating a pool of potential trachea candidates, TCP1, in VOI1; and on condition that there is one 3-dimensional, 3D, component in the TCP1, locating the trachea as the 3D component; on condition that there is more than one 3D component in the TCP1, determining a second volume of interest, VOI2; segmenting air areas in VOI2; generating a second pool of potential trachea candidates, TCP2, in VOI2; combining TCP1 and TCP2 to find common 3D components; locating the trachea as the 3D component spanning the largest total number of CT images of the combined TCP1 and TCP2; and selecting a slice from the CT volume of images where the trachea is located as a seed point S_t for a subsequent extraction of the trachea, wherein detecting (S12b) any tracheal intubation present comprises: calculating a center point of the trachea in a binary image from the image stack; performing image thresholding on the binary image; excluding areas in the binary image having less than 15 and above 100 pixels; applying 2D connected component analysis to the binary image; calculating center points of all distinct 2D connected components in the binary image; calculating a Euclidean distance of component center points to the center point of the trachea to identify the closest components to the trachea center point; identifying any center point of the 2D connected components having a distance of less than 11 pixels as a potential endotracheal tube; and on condition that a tracheal intubation is present, performing 3D region growing to the CT data starting at the seed point S_t to extract the tracheal intubation; and identifying a slice in the image stack of the extracted tracheal intubation in cranial to caudal direction as the last image where the tracheal intubation exists, and selecting said slice as seed point S_t for a subsequent extraction of the trachea, wherein segmenting (S12c) the airways comprises: thresholding the selected slice to obtain a binary image; and segmenting and tracking the trachea region up to bifurcation point in cranial to caudal direction using the 2-dimensional, 2D, growing method using the binary image and the seed point S_t as inputs, wherein the selected slice and seed point S_t is taken from step 12a if no tracheal intubation exists, and from 12b if tracheal intubation exists, wherein identifying (S12d) the optimal carina location as the carina trachea comprises: identifying the optimal carina location by tracking the trachea slice by slice from the bifurcation point in cranial to caudal direction, calculating the distance between left main bronchus and right main bronchus for each respective slice until reaching a slice, Soc, where the calculated distance is equal to or greater than 0,75 cm, and selecting said slice Socas the optimal carina location; and identifying the selected optimal carina location as the carina trachea. The method of any one of claims 1-7, wherein detecting (S13) descending aorta at carina level comprises: creating a first search space originating from an artificial line that cuts the x- axis of the slice at an angle of 22,5 -75 ° downwards clockwise from the center of left main bronchus comprising lines having an angle of angle larger than 2,5° and smaller than 15° between each other forming rays from center of origin; creating a second search space originating from an artificial line that cuts the x-axis of the slice at an angle of 60 -120 ° downwards clockwise from the center of the left main bronchus comprising lines having an angle larger than 2,5° and smaller than 15° between each other forming rays from center of origin; identifying if any one of the rays pass through the descending aorta by going ray by ray and: collecting all pixels on the ray; thresholding said pixels over - 1 HU to acquire binary image; applying 2D connected component analysis to identify largest connected pixels on the ray; calculating the total number of pixels on the segmented component; and if the total number of pixels is greater than 15; calculating mean x and y coordinates of those pixels, separately, in order to reduce to one data point, and identifying said data point as the descending aorta at carina level, PdesAo. The method of any one of claims 1-8, wherein detecting (S14) ascending aorta comprises: detecting (S14a) aortic arch by tracking the descending aorta up to the aortic arch in caudal to cranial direction; detecting (S14b) the ascending aorta by tracking the upper part of the aortic arch to the ascending aorta in cranial to caudal direction to locate a first slice of the ascending aorta. The method of any one of claims 1-9, wherein detecting (S16) pulmonary trunk comprises: creating (S16a) a rectangular search space adjacent to the left lateral side of the ascending aorta location in the CT image stack covering the candidate of the pulmonary trunk, the pixel location of the rectangular being calculated as (eAot - 15) eAob, (E4OZ +5) : eAoi +70)), where eAo_t is the top most and eAob the bottom most pixel of the ascending aorta region; and detecting (S16b) the pulmonary trunk within said rectangular search space. The method of any one of claims 1-10, wherein detecting (S17) an apical level of pulmonary valve comprises: segmenting and tracking (S17a) the pulmonary trunk slice by slice from carina level to conus arteriosus in cranial to caudal direction; identifying (S17b) a region segmented in the last slices from step 17a; evaluating (S17c) sphericity of the segmented region using a first criterion, ciri, to identify the apical level of pulmonary valve; on condition that the apical level pulmonary valve is not identified using the first criterion, evaluating (S17d) the sphericity of the segmented region using a second criterion, c+2, to identify the apical level of pulmonary valve, detecting (S17e) an apical level of the pulmonary valve as the last segmented component identified in step S16b. A computer-aided detection, CADe, method performed in a computing device for locating and measuring large vessels and detecting ascending aorta hypertension and/or pulmonary trunk hypertension in a subject, the method comprising: performing the method of any one of claims 2-11 to determine a representative diameter of an ascending aorta and/or pulmonary trunk in a subject; comparing the determined representative diameter(s) with a preset threshold; determining that ascending aorta hypertension and/or pulmonary trunk hypertension is present in a subject if the determined respective representative diameter(s) is greater than the preset threshold. The method of claim 12, wherein a presence of ascending aorta hypertension and/or pulmonary trunk hypertension in the subject indicates presence of a medical disorder in the subject, the medical disorder being one or more of pulmonary embolism, heart disease, thrombo-embolism, cryptogenic organizing pneumonia and scleroderma. A computer-aided detection system (200) comprising a computing device (10), configured for performing a computer-aided detection, CADe, method of locating and measuring large vessels in a subject undergoing a computed tomography, CT, scan, the computing device (10) comprising: a memory (13) for storing instructions and processing circuitry (12) configured to cause the computing device (10): to obtain CT data (20) from a consecutive CT scan of thorax of the subject, the CT data (20) comprising an image stack forming a CT volume of images; to prepare the CT data (20) for segmentation tasks; to locate carina trachea in the CT volume; to detect descending aorta at carina level; to detect ascending aorta; to detect pulmonary trunk; and to detect an apical level of pulmonary valve. The computer-aided detection system of claim 14, configured to perform the method of any of claims 1-11. A computer program comprising computer program code which, when executed in a computing device, causes the computing device to execute the methods according to any of the claims 1-11.