Optical Flow-Based Extraction of Breathing Signal from Cone Beam CT Projections

Abstract

Respiratory motion serves as a major challenge during treatment of lung cancer patients using radiotherapy. In this work, an image-based method is presented to extract a respiratory signal directly from Cone Beam CT (CBCT) projections. A dense optical-flow method is used to acquire motion vectors between successive projections in each dataset, followed by the extraction of the dominant motion pattern by application of linear kernel Principal Component Analysis (PCA). The effectiveness of the method was tested on three patient datasets and the extracted breathing signal was compared to a ground-truth signal. The average phase shift was observed to be 1.936 ± 0.734 for patient 1, 1.185 ± 0.781 for patient 2 and 1.537 ± 0.93 for patient 3. Moreover, a 4D CBCT image was reconstructed, considering the respiratory signal extracted, using the proposed method, and compared to that reconstructed considering the ground-truth respiratory signal. Results showed that a minimal difference was found between the image reconstructed using the proposed method and the ground-truth in terms of clarity, motion artifacts and edge sharpness.

1. Introduction

Lung cancer cases are on the rise globally, with 2.21 million reported cases in the year 2020. It is also a top contributor to cancer-associated mortality as per the statistics provided by World Health Organization [1]. Lung cancer treatment involves radiotherapy in which a radiation pulse is targeted at the tumor site to destroy cancerous cells. The treatment planning process during radiotherapy for cancers of the lung and superior abdominal regions is affected by the breathing motion. Accounting for the breathing motion is necessary through motion management techniques to ensure that the location of the tumor is not misrepresented, allowing the radiation beam to precisely kill the cancerous cells and at the same time minimize the exposure of healthy tissues surrounding the tumor site to a high dose of radiation.

Several motion management techniques for application in radiotherapy have been developed and tested [2]. Among such techniques commonly used are the application of surface markers that track chest motion as it moves during respiration or the implantation of internal fiducial markers that can help track the tumor motion with better accuracy [3,4,5,6]. It is not evident how well a signal estimated by tracking the movement of a marker block positioned on the skin correlates with a tumor located within the lung. On the other hand, implantation of fiducial markers can pose health risks to the patient. Use of techniques such as breath hold and abdominal compression can minimize the effect of motion on the images; however, such techniques can cause discomfort to the patient undergoing the scan.

Therefore, there is a need to shift towards image-based methods for management of breathing motion that does not rely on marker tracking, and are robust and safe. Image-Guided Radiation Therapy (IGRT) is an emerging technology adopted in clinical practice due to its ability to deliver highly accurate radiation doses to the tumor and reduce the side effects at the same time. Based on a master’s thesis [7], this paper presents an image-based method to estimate respiratory motion using Cone Beam Computed Tomography (CBCT) projections. In this work, an optical flow-based approach is first used to produce motion vectors. This process is followed by the application of Principal Component Analysis (PCA) to extract dominant motion patterns from the flow vectors. Moreover, the extracted respiratory signal is used to sort the CBCT projections according to their respiratory phase, which can be used to reconstruct a 4D CBCT image. The reconstructed image is evaluated by comparing it to the image reconstructed considering the ground-truth breathing signal. The contributions of this work lie in building an image-based method that is simple, computationally efficient, and can be used for respiratory phase sorting and 4D CBCT image reconstruction right after the CBCT projections are acquired.

The subsequent sections of this paper are organized in the following manner. Section 2 includes the literature review. Section 3 discusses the used datasets and methods. Section 4 presents and discusses the experimental results. Section 5 provides a discussion of the results and limitations of this work, and provides future research directions. The paper is concluded in Section 6.

2. Literature Review

Several studies are found within the literature that show the potential of image-based techniques, and some of the relevant studies are briefly discussed in this section. Figure 1 shows some of the commonly used respiratory motion-management techniques, including the image-based methods.

2.1. Slow CT Scanning

Slow CT scanning is a method to encompass tumor volume. During the slow CT scanning procedure, the scanner rotates at a low speed of about 4 s per slice, such as to match the duration of an average respiratory cycle. Averaging of several CT scans can also be performed to capture multiple breathing phases in each slice. The phantom study performed in [8] demonstrate that with slowing down the period of gantry rotation to 1.5 s, it is possible to better reproduce and identify target motion. During treatment planning, a slow CT scan aids in outlining tighter margins for the target volume [9]. This method is generally not advisable for tumors located near the diaphragm or chest wall, because of the loss in resolution that can occur due to blurring resulting from motion.

2.2. Breath Hold

During the breath-hold technique, the patient holds their breath for specific intervals, during which images are acquired. In comparison to scans taken with slow CT, images under breath-hold are less blurry, with reduced artifacts. M. J. Murphy et al. [10] performed a study to evaluate the effectiveness of the breath-hold technique in replicating the tumor location during radiotherapy. It was possible to reproduce the tumor location at an average of within 1 millimeter for the thirteen patients included in the study. It was also highlighted that the motion margin for dose delivery can be lowered to less than 2 millimeters with the breath-holding technique. In order to assist patients with breath-hold, active breath coordinator (ABC) devices can be used. The most prominent limitation of this method is that it causes discomfort to the patient, as it requires them to hold their breath for at least 10 s.

2.3. Abdominal Compression

Respiratory motion management through the abdominal compression technique involves the application of pressure on the abdomen through compression plates. J. H. Heinzerling et al. [11] tested the effect of different degrees of abdominal compression on the movement of tumor in the lung. For no compression, medium compression and high compression, the average tumor motion was observed to be 13.6 mm, 8.3 mm and 7.2 mm, respectively. Although the overall movement can be reduced by applying abdominal pressure, this technique causes inconvenience to the patient undergoing the scan.

2.4. Surface Block Markers

The breathing signal during radiotherapy can also be estimated through external surrogates such as block markers and spirometer. The blocks contain infrared markers and are placed on the abdomen or chest of the patient. The movement of the block in superior–inferior direction is captured by infrared camera. H. Yan et al. [12] performed a study to evaluate how the correlation between a block marker placed externally on the patient’s body and the actual tumor motion is affected by the placement location of the marker. For the seven patients included in the study, the correlation coefficient ranged between 0.247 and 0.987. The research highlighted the fact that there can be several factors that can influence this correlation, such as characteristics of the patient, pattern of breathing and position of the marker on the skin. The marker position can be patient-specific; therefore, it needs to be adjusted, based on the individual. The signals from multiple markers are often combined to better track the tumor motion.

2.5. Spirometer

A spirometer can also serve as an external surrogate for tumor motion. The spirometer represents the lung tidal volume, providing a measure of the amount of air inhaled and exhaled. The change in lung volume assessed with the spirometer can be correlated with tumor location [13]. Spirometer-based estimation of the respiratory signal provides an advantage over a respiratory surrogate obtained by tracking of surface block markers. Placement of the infrared blocks and camera angle can vary between different sessions, which can give rise to error while estimating the tumor position. However, with spirometer procedures there is no change in the everyday setup. The disadvantage of this method is the baseline signal drift, which makes it inefficient for use in long treatment procedures.

2.6. Image-Based Methods

2.6.1. Fiducial Marker Tracking

Fiducials are in the form of gold seed markers placed surgically near the tumor site, and are about 2 to 4 mm in diameter. Fiducials help in the better localization of the radiation beam to the target area, and protect the healthy tissues from harmful radiation exposure. One of the essential requirements for accurate monitoring of tumor motion through fiducial tracking is that the fiducial does not migrate from its implanted location over time. M. Imura et al. [14] investigated fixation rate and displacement of a 1.5-millimeter marker implanted in lung cancer patients. About 75% of the fiducials were identified for the entire treatment process. It was observed that the markers are better attached in bronchial lumens having narrower diameters. The internal marker and tumor attachment suffered significant variation after a period of two weeks. Implantation of metal markers poses certain risks for the patients, such as causing the lung to collapse, known as pneumothorax, spitting of blood (hemoptysis) or haemorrhage [4]. N. Kothary et al. [15] studied the safety aspect of a fiducial marker used during radiotherapy. Out of the 44 markers implanted in lungs of different patients, 20 were observed to have shown pneumothorax. Even though internal markers are effective in localizing the tumor site during radiotherapy, the associated health risks raise the need to shift towards markerless estimation of respiratory motion.

2.6.2. Diaphragm Tracking

Image-based extraction of the respiratory signal can also rely on the diaphragm as an anatomical marker. The diaphragm is visible as a region of high contrast in the scanned images, distinguishing the lungs from the upper abdominal region. In the study by L. I. Cervino et al. [16], the correlation between diaphragm movement and motion of the tumor using two linear models was investigated. The main difference between the two models was that the second model took into consideration the shift in the phase between tumor and diaphragm movements. A mean correlation coefficient of 0.94 and 0.98 were observed with the first- and second-motion models. Even though a good correlation factor was observed with both models, it was emphasized that the correlation can differ between different patients and it must be investigated for each patient individually. J. Wei and M. Chao [17] performed a study to extract respiratory motion from CBCT images through diaphragm tracking. A medical professional determined five points on the projection images to identify the area containing the diaphragm. Through a modified canny edge detection scheme, the identified points were detected in the successive images throughout the projection data. In comparison to manual identification of the diaphragm edge, the observed error with the proposed algorithm for CBCT projections was 0.79 mm. For the rotating cone-beam CT, the diaphragm shape differs at different angle of gantry rotation, making it difficult to reliably track diaphragm position across the entire set of projections.

2.6.3. Markerless (Amplitude/Phase) Binning

Respiratory motion-management methods such as surface markers, abdominal compression, breath-hold, fiducial marker and diaphragm tracking face certain limitations. While it is not very evident how well a marker placed on patient’s chest correlates with a tumor motion within the lung, fiducial markers implanted near the tumor site can cause certain health risks. Breath-hold and abdominal-compression methods try to minimize the motion, but cause discomfort to the patient. To overcome the limitations caused by these methods, a shift towards markerless image-based methods is needed, to estimate the breathing motion and reconstruct motion-compensated 4-dimensional images. Several studies are found within the literature that have highlighted the potential of using markerless methods.

Y. Kubota et al. [18] designed an algorithm using an improved 3D KLT (Kanade–Lucas–Tomasi) tracker to estimate lung motion from 4D CT scans. The tracker identified feature points from the projections based on change in image gradient or intensity. It was observed that the feature points near the inferior area of the lung show more displacement and provide better tracking accuracy as compared to feature points in the superior lung region. The tracking accuracy is also significantly reduced with higher rotation angle, going from an accuracy of 79.6% at 2 degrees to 15.8% at 10 degrees.

I. Vergalasova et al. [19] investigated the feasibility of a novel markerless technique to sort CBCT projections for 4D reconstruction. As per the Fourier transform (FT) principle, the FT image magnitude shows how much a certain frequency component is exhibited in an image; on the other hand, the phase shows the position of that specific frequency in the image. Projections associated with the peak inhale phase were identified with both amplitude and phase FT methods. The methods were tested on two phantom and three clinical datasets. For the FT phase method, the average phase shift compared to manual ground-truth for phantom 1 with 3 s breathing cycle was 1.8%, and for phantom 2 with 6 s breathing cycle it was 3.9%. The clinical datasets resulted in a difference of 2.9%, 5.0% and 3.8% for the three patients, respectively. The magnitude FT method resulted in a difference of 2.1% and 4.0% for the two phantoms and 2.9%, 5.3% and 3.5% for the three patients. The observed phase difference and quality reconstructed images using both methods show the feasibility of Fourier transform methods. More investigation on patient datasets attained at slow gantry rotation is required to be performed, to fully understand the potential of the proposed technique.

S. Dhou et al. [20] developed a markerless technique for the extraction of the respiratory signal called local intensity feature tracking (LIFT). Firstly, feature points were extracted at evenly distanced pixels, excluding the pixels representing the diaphragm region. The feature points were then tracked across the projection datasets and motion trajectory was formed. Trajectories experiencing oscillatory behavior similar to the respiratory pattern were clustered, thereby retrieving the 3D lung movement. A mean phase shift of 1.68 projections and 1.78 projections was observed with LIFT, in comparison to ground-truth obtained from diaphragm tracking and internal marker tracking, respectively. S. Dhou et al. [21] suggested the use of the optical-flow technique with linear and non-linear dimensionality reduction to extract the respiratory signal. Implementation of the method on patient data resulted in a mean phase shift of 3.77 ± 1.31, 3.59 ± 2.44, 3.83 ± 2.65 for the three patients 1, 2 and 3, respectively.

S. Park et al. [22] presented an image registration-based projection binning (IRPB) method for sorting 4D CBCT projections. The proposed method involved identification of feature points in accordance with changes in intensity inside a rounded pixel mask. Trajectories of feature points over the sequence of CBCT projections were then formed, using a random sample consensus-image-registration method, followed by the application of PCA, to identify the trajectory related to the respiratory signal. In comparison with other image-based methods for breathing-signal extraction, such as the Fourier transform method and LIFT, the average phase shift with this technique was observed to be 3.74 and 0.48 projections lower, respectively. Peak signal-to-noise ratio for reconstructed images was 5.08 dB, 1.05 dB and 2.90 dB higher than with the Amsterdam Shroud, Fourier transform and LIFT methods, respectively. cIRPB method provide several advantages, such as being markerless, and overcomes the restrictions caused by limited-view projections. The authors highlighted the fact that their method is limited by the processing time and memory requirements.

Another promising approach for extracting lung tumor motion from CBCT projections was presented by M. Chao et al. [23]. Initially, planned CT projections were separated into two groups, with the first group of projections having only the tumor and the second group having the tumor removed. Rigid frame-wise registration was then carried out between CBCT projections and digitally reconstructed radiographs (DRRs) of the two groups. Subtraction of the registered DRR from CBCT projections resulted in intensified tumor visibility. Tumor motion was tracked through active contouring. In comparison with the ground-truth signal, the signal extracted through the proposed method resulted in an error of lower than 0.2 mm. The error for patient data was between 1.4 and 2.2 mm.

In addition to the mentioned studies, other significant studies in the literature include P. Fischer et al. [24], who used PCA on fluoroscopy images to extract breathing signal. T. Geimer et al. [25] used kernel-ridge regression on planning and follow-up 4D CT data. In this work, we utilize the CBCT data only, as it provides the most recent view of the internal anatomy. A preliminary version of our research performed on phantom datasets has been previously reported and published [26].

The emerging artificial intelligence (AI) techniques have recently been adopted in different domains. AI has been applied in IGRT and has shown great potential in motion tracking in radiotherapy. In a recently published review paper [27], four main categories of motion management using AI are examined and summarized: (1) marker-based and (2) markerless tracking, where, in these two categories, the individual target is tracked throughout the treatment, (3) full-anatomy monitoring, where the intrafraction changes are monitored in the full anatomy within the field of view, and (4) motion prediction, where AI is used to account for the time latencies to localize, process and act.

3. Materials and Methods

3.1. Datasets

The present study is performed on patient CBCT datasets, as an extension to the study that was previously performed on phantom datasets [26]. Three anonymous patient datasets were included in the study. The patient datasets were shared by Virginia Commonwealth University, Richmond, Virginia, USA, under a Data Transfer and Use Agreement, and were retrospectively analyzed. Moreover, the research protocol qualified for exempt approval from the Institutional Review Board (IRB) of the American University of Sharjah, United Arab Emirates (Approved IRB Protocol #: 18-425). The patient datasets consist of CBCT scans that were performed by Elekta and Varian CBCT imaging systems. The specifications of the patients datasets can be found in

Table 1. Description of the patient datasets.

Dataset	Projection Size	Number of Projections	Imaging System	Ground-Truth Breathing Signal
Patient 1	512 × 512	701	XVI 3.5 (Elekta)	Diaphragm motion
Patient 2	768 × 1024	1220	OBI (Varian)	Internal marker trace
Patient 3	768 × 1024	2396	OBI (Varian)	Internal marker trace

.

3.2. Feature Tracking Using Optical Flow

Optical flow is the distribution of 2 D velocities in an image caused by the motion of luminance patterns also known as pixel intensities [28]. The dense optical-flow method with enhanced classical Horn–Schunck method [28] provided by D. Sun, S. Roth, and M. J. Black [29] is implemented in this work. For two subsequent frames acquired at time $t$ and $t+1$ in a sequence of images, pixel $\left(x,y\right)$ in frame 1 changes position by $\left(x+\delta x,y+\delta y\right)$ in frame 2. The motion of pixel intensities in an image sequence can thus be described by Equation (1).

$$I\left(x+\delta x,y+\delta y,t+1\right)=I\left(x,y,t\right)$$

(1)

The apparent motion of pixel intensities between subsequent projection images while estimating respiratory motion can be given as Equation (2).

$${\displaystyle \frac{\partial I}{\partial x}}{\displaystyle \frac{dx}{dt}}+{\displaystyle \frac{\partial I}{\partial y}}{\displaystyle \frac{dy}{dt}}+{\displaystyle \frac{\partial I}{\partial t}}=0$$

(2)

Let $U={\displaystyle \frac{dx}{dt}}$ and $V={\displaystyle \frac{dy}{dt}}$ be the horizontal and vertical optical-flow components of I(x,y,t).

Then,

$$I_{x}U+I_{y}V+I_t=0$$

(3)

where $I_x, I_y, I_t$ represent the partial derivative of image intensity in respective directions. Motion fields were estimated through a coarse-to-fine approach. For every adjacent projection image, an image pyramid is built by sub-sampling. The base level of the pyramid consists of the original frame of size $\left(M\times M\right),$ and adjacent layers are occupied as $(1\times 1, 2\times 2, 4\times 4,\dots , 2^a\times 2^a)$ so that $M=2^a$. The pyramidal approach involves the computation of optical flow initially at the layer with images at minimum resolution, to the layer with the original frame.

In this work, the optical flow is applied to all image pixels of every projection. Thus, a data matrix is created for each dataset to store the output of the optical-flow vectors. The total number of columns in the matrix is equal to (2 × the dimension of the projection), as each pixel is described by two components: the horizontal ${∆u}_{\left(prj 1, pxl 1\right)}$ and the vertical one ${∆v}_{\left(prj 1, pxl 1\right)}$. On the other hand, the number of rows in the matrix is equal to the total number of projections. The matrix is filled row by row, such that all pixel displacement components for the first CBCT projection are filled in the first row of the matrix, and this is continued till the final CBCT projection.

The data matrix $D$ can be described by Equation (4).

$$D=\left[\begin{array}{c}{∆u}_{\left(prj 1, pxl 1\right)}\\ ⋮\\ {∆u}_{\left(prj n, pxl 1\right)}\end{array} \begin{array}{c}{∆v}_{\left(prj 1, pxl 1\right)}\\ ⋮\\ {∆v}_{\left(prj n, pxl 1\right)}\end{array} \begin{array}{c}\cdots \\ ⋮\\ \cdots \end{array} \begin{array}{c}{∆u}_{\left(prj 1, pxl m\right)}\\ ⋮\\ {∆u}_{\left(prj n, pxl m\right)}\end{array}\begin{array}{c}{ ∆v}_{\left(prj 1, pxl m\right)}\\ ⋮\\ {∆v}_{\left(prj n, pxl m\right)}\end{array}\right]$$

(4)

where prj and pxl represent projection number and pixels, respectively. The final projection within a set of data is represented as n, and m represents the last pixel in the projection being processed. $∆u$ and $∆v$ signify the horizontal and vertical component of the optical-flow vectors.

3.3. Respiratory Motion Analysis Using Principal Component Analysis (PCA)

In the next stage, respiratory motion is extracted from the optical-flow data matrix by application of Principal Component Analysis (PCA) with linear kernel. PCA is an unsupervised learning technique which mathematically translates interrelated variables to principal components by projecting data into another coordinate system. PCA is considered in this work due to its efficiency, simplicity and ability to distill the large amount of motion vectors into a set of a few eigen vectors representing the breathing patterns in the dataset.

In this work, PCA is performed with the help of a dimensionality reduction toolbox [30]. The process involves computation of a gram matrix, as given by Equation (5).

$$K=D\ast D^{\prime }$$

(5)

Next, eigenvalues ($\lambda )$ and eigenvectors ($V)$ of the normalized K matrix are computed, and are sorted in decreasing order. The output of PCA can be described by Equation (6):

$$L=[{\lambda }_1 {\lambda }_2 \dots {\lambda }_p] , T=\left[\begin{array}{ccc}{vect}_{\mathrm{1,1}}& \cdots & {vect}_{1,p}\\ ⋮& ⋮& ⋮\\ {vect}_{prj,1}& \cdots & {vect}_{prj,p}\end{array}\right]$$

(6)

where p represents the number of the principal component, T represents the eigenvector matrix and L represents eigenvalues for every principal component along the diagonal. The signal from every principal component is computed as $\left(V^{\prime }\ast diagonal of \sqrt{\lambda }\right)$. The principal eigenvectors represent the prominent or main variations in the flow vectors, and hence the motion of lungs during breathing can be extracted by linear combination of eigenvectors exhibiting maximum eigenvalues.

3.4. Respiratory Phase Sorting and 4D CBCT Reconstruction

The next step after acquiring the respiratory signal involves phase sorting, which aims at binning the projections into the corresponding respiratory phase. The phase-sorting process is presented in Figure 2. For each breathing cycle, phase 1 is assigned to the max exhale. The overall number of projections in the breathing signal is divided by the total number of phases, and the remainder of the projections are assigned phase numbers, accordingly. The division is performed such that the count of projections for each phase bin is roughly similar.

The accuracy of the extracted breathing signal and the phase-sorted signal is tested by performing a comparison with ground-truth and calculating the average phase-shift error. The ground-truth for the datasets includes the diaphragm position and internal marker traces.

The main of aim of respiratory motion extraction and phase-sorting procedures is to obtain projections corresponding to the same phase, which can be used to reconstruct 4D CBCT images. Projections belonging to the similar phase exhibit very little or no motion. Four-dimensional reconstruction in this work was performed through the FDK algorithm with the help of the RTK reconstruction toolkit [31]. The reconstructed 4D images can be assessed by plotting the edge profiles to evaluate the sharpness of organ edges in reconstructed slices, and the calculation of the reconstruction error and peak signal-to-noise ratio (PSNR).

4. Results and Discussion

In this section, the experimental results are presented. The MathWorks, Inc., Natick, MA, USA is utilized as the main programming software in this research.

Results for patient 1 dataset:

This dataset comprises 701 CBCT projections with ground-truth obtained by diaphragm tracking. Figure 3 shows the extracted breathing signal for the patient 1 dataset and the corresponding phase-sorting plot. The average phase shift and standard deviation was observed to be 1.936 and 0.734, respectively.

Results for Patient 2 dataset:

This dataset comprises 1220 CBCT projections of a dimension of 768 × 1024 with ground-truth obtained by tracking the location of one implanted marker. Figure 4 shows the extracted breathing signal for two different projections from the patient 2 dataset and the corresponding phase-sorting plot. As the total number of projections was too large to be processed at one time, the processing was therefore performed for groups of 300 projections. The average phase shift and standard deviation phase shift of the entire signal for the patient 2 dataset was observed to be 1.185 and 0.781.

Results for Patient 3 dataset:

This dataset comprises 2396 CBCT projections with a dimension of 768 × 1024, with ground-truth obtained by tracking the location of four internal markers. Similar to the patient 2 dataset, the processing was performed for a group of 300 projections at a time. A plot of the extracted breathing signal in red, along with the position of four markers and their average location, is shown in Figure 5. The extracted signal was able to capture sudden changes in respiratory pattern similar to what was captured by markers A and B. These changes are highlighted in the plot with a black dotted rectangle. Figure 6 shows the extracted breathing signal for two different projetions from the patient 3 dataset and the corresponding phase-sorting plot. The average phase shift and standard deviation phase shift of the entire signal for the patient 3 dataset was observed to be 1.537 and 0.93. Unlike signals captured by markers C and D, the signal extracted by the proposed algorithm is fully continuous.

The average phase shift for patient 1 is higher in comparison to patient 2 and patient 3. The main difference between these datasets is that for dataset 2 and 3 the ground-truth was provided through an internal marker, while for patient 1 it was through diaphragm tracking. This shows that the extracted signal correlates more with the extracted signal using the internal markers, as they are both extracted from internal lung tissues.

The extracted breathing signal was used for phase sorting where projections belonging to the same breathing phase were used for 4D image reconstruction. To test the feasibility of using the extracted signal for phase sorting and 4D image reconstruction, the proposed breathing signal extraction method was applied to a phantom dataset. This dataset was generated using the XCAT (extended cardiac torso) simulation software by W. P. Segars et al. [32,33]. The dataset was prepared simulating regular respiratory motion, where every cycle comprises 20 breathing phases. The total number of projections in this dataset is 720. The respiratory motion extraction results for this dataset were previously published in [26], where the dataset was referred to as Phantom dataset 2. Figure 7 shows coronal slices from the 3D and 4D CBCT reconstructed images of the phantom dataset.

As is shown in Figure 7, when all phases are used for reconstruction, as in (a), the organ edges appear blurry, while they become significantly more distinct when reconstructed with a single phase obtained after binning, as in (b) and (c). In this figure, phase 1 was considered. In (b), the projections were sorted considering the ground-truth breathing signal, while in (c) the projections were sorted considering the extracted breathing signal, using the proposed method. The image in (d) is the result of subtracting the images in (b) and (c). As can be seen in (d), the difference is mainly in demonstrating noise and streaking artifacts in the image, and is unrelated to structural or anatomical differences in the images.

To further evaluate the reconstructed images, edge profiles extracted from each of the reconstructed images were compared. The edge profile helps to evaluate the sharpness of organ edges in reconstructed images. The profile considered is extracted from the same slice and from the same location in each of the reconstructed images. It can be demonstrated by a yellow vertical line, shown in Figure 7a. Figure 8 shows the edge profiles extracted from each of the reconstructed images, namely the 3D reconstructed image using all projections, as shown in Figure 7a, the one reconstructed from phase 1, which only considers the ground-truth breathing signal, as shown in Figure 7b, and the one reconstructed from phase 1, which only considers the breathing signal extracted using the proposed method, as shown in Figure 7c. The results in Figure 8 are in accordance with the visual results provided in Figure 7, where the organ edges appear blurry when all phases are used for reconstruction, while they become significantly sharper when reconstructed with a single phase. Furthermore, the mean square error for the image reconstructed considering the breathing signal extracted using the proposed method was observed to be 2.65 × 10⁻⁵, while the PSNR was observed to be 45.75.

5. Discussion

The experimental results showed that the proposed method was able to extract the respiratory motion from the CBCT images with a minimal difference compared to the ground-truth signal. The work presented in this paper can be compared to the works presented in [20,21,26]. In [20], the average phase shift was reported as 1.68 projections when compared to signal extracted based on diaphragm position, and 1.78 projections when compared to signal extracted based on internal markers. The work presented in [20] involved several steps which may increase the method complexity. In our previous work in [26], the average phase shift did not exceed 1.5 projections for tests performed on regular and irregular phantom datasets. However, that method was only applied to phantom datasets. In [21], the average phase shift was 3.8 projections when compared to signal extracted based on diaphragm position, and 3.59 projections when compared to signal extracted based on internal markers using the non-linear dimensionality reduction method which was the best performing method. Although applied to clinical datasets, the average phase shift was relatively high. Furthermore, no 4D CBCT reconstruction results were provided. All in all, compared to [20,21,26], the extracted respiratory signal using the proposed method in this paper achieved lower phase shift than similar results in the literature.

Moreover, the reconstructed 4D CBCT images considering the respiratory signal extracted using the proposed method showed clear images with sharp edges compared to the ground-truth based images. The reconstructed 4D CBCT images are affected by streaking artifacts that are caused by the reconstruction process wherein only un-der-sampled projections are used for the reconstruction of each phase. These artifacts have nothing to do with the phase sorting process. Several studies have been proposed in the literature to minimize the streaking artifacts and improve the overall quality of the reconstructed 4D CBCT images. These solutions revolve around using motion compensation and modeling methods [34,35], deep learning [36,37,38] or projection oversampling through in-between projection interpolation [39,40].

The limitations of this work include the limited number of datasets used to test the proposed method. As a future work, it is recommended to test the proposed method on a larger number of datasets to get more reliable findings. Moreover, different optical flow methods can be tested for motion estimation. In addition, other dimensionality reduction methods can be explored for the breathing motion extraction from the optical flow vectors

6. Conclusions

In this work, an algorithm for the extraction of a breathing signal for the purpose of phase binning and 4D CBCT reconstruction was provided. The algorithm is image-based, and does not rely on any sort of surface marker or implanted fiducials. Experiments were performed on several patient datasets The experimental results showed signals correlated well with the ground-truth signal, with the highest average phase shift of 1.936 ± 0.734 being observed for patient 1 when the ground-truth was provided through diaphragm location and not the internal marker. Four-dimensional CBCT reconstruction was also performed on the phantom dataset to show the effectiveness of the proposed algorithm. The present work can be further expanded to include more patient data in the study and test other optical flow and dimensionality reduction algorithms to perform a comparison.