Remote monitoring of breathing dynamics using infrared thermography

Carina Barbosa Pereira; Xinchi Yu; Michael Czaplik; Rolf Rossaint; Vladimir Blazek; Steffen Leonhardt

doi:10.1364/BOE.6.004378

1. Introduction

Breathing rate (BR) is considered an essential vital sign [1, 2]. In adults, this parameter varies from 12–18 bpm (breaths per minute) under resting conditions [3], whereas neonates’ BR is ranging between 30 and 50 bpm [4]. An abnormal breathing manifests itself as tachypnea (high BR), bradypnea (low BR), phases of apnea (transient cessation of breathing) or as irregularities in breathing. Moreover, some pathologies (e.g. metabolic disorders [5] and cerebral ischemia [6, 7]) lead to specific breathing patterns such as Kussmaul breathing and Cheyne-Stokes respiration, respectively. Among others, occurring breathing pathologies serve as solid and one of the earliest indicators of physiological distress [1, 2]. In addition to that, monitoring of breathing plays a crucial role in diagnosing and managing breathing disorders such as sleep obstructive apnea, bronchitis, asthma, etc. Sport and sleep studies, as well as early detection of sudden infant death syndrome (SIDS), are further examples where BR is thought to be a highly relevant parameter [3]. As a matter of fact, SIDS is still one of the major causes of death in infants younger than 1 year.

Nevertheless, BR is often a neglected and underestimated parameter mostly due to shortcomings of the established monitoring techniques [1, 2, 8, 9]. They require attachments of sensors to the patient, leading to discomfort, stress and even to soreness, particularly in small children [10]. In general, monitoring of BR can be performed by (1) measuring chest or chest-abdomen movement with respiratory belt transducers [11, 12], (2) detecting acoustic signals of the upper airway which occur due to turbulent airflow [13–15], (3) metering electrical impedances of the thorax (e.g. by using electrical impedance tomography - EIT) [16, 17] or(4) performing signal analysis from electrocardiography (ECG) [11, 18, 19]. Moreover, there are more precise but also more invasive methods available such as capnography [19], monitoring varying partial pressure of carbon dioxide (CO₂) in tidal volume, as well as spirometers [11] and nasal thermistors [3] measuring air-flow and nasal temperature variations, respectively.

In recent years, there has been an increasing demand for unobtrusive and contactless but also reliable monitoring alternatives of BR, aiming at improving patients’ quality of life and optimizing the use of medical resources [1]. Therefore, new monitoring solutions based on Doppler radar [1] and imaging sensors (visible [11], midwave infrared [12] and long wave infrared imaging sensors [10, 20]) have been being proposed and developed.

Thermal imaging, also denominated infrared thermography (IRT), emerged as a promising monitoring and diagnostic technology in a wide range of medical fields. Heart rate [21] and BR detection [10, 12, 22], monitoring of thermoregulation in neonates [23] and observation of circulation and perfusion dynamics [24] are some examples where the application and performance of IRT was studied. Thermal imaging is a remote non-contact monitoring approach. In addition, it is a passive technique, i.e. it detects the radiation naturally emitted from an object, in this case the human skin, and does not use any harmful radiation. Last but not least, thermal imaging does not need a light source. This particular characteristic is one of the biggest advantages of IRT over other imaging technologies [24].

This scientific paper proposes a novel reliable approach for contactless and passive monitoring of breathing function in thermal imaging. In contrast with other methods, which need a manual selection of the region of interest [12, 22, 25], our approach allows an automatic detection of the nose in the first frame. In addition, a flexible algorithm for robust estimation of local breath-to-breath intervals from IRT waveform was implemented. It may be a good alternative to extract BR in nonstationary signals. While other research groups just validated the algorithm in optimal conditions [12, 22] (i.e. the participants were instructed to minimize head movements and to breathe normally), in this study we address the ability of our approach to accurately measure breathing rate during challenging conditions: head motion and breathing disorders. Here, its robustness against motion artifacts as well as its accuracy and reliability during breathing maneuvers were tested. The developed methodology, particularly segmentation and tracking of the nostrils, extraction of the breathing waveform as well as its processing, is described in section 2. Section 3 introduces, in turn, the experimental setup. Results are presented and discussed in section 4 and section 5, respectively. Finally, section 6 concludes the paper and gives future perspectives.

2. Methodology

The current approach is based on the fact that temperature around the nostrils fluctuates during the respiratory cycle (inspiration and expiration). Whereas during inhalation cold air from the environment is inspired, during expiration warm air from the lungs is exhaled. IRT is capable of accurately detecting this nasal temperature modulation as depicted in Fig. 1. In this example, the mean temperature of the region of measurement (ROM) corresponds to 31.17 °C during inspiration and 31.44 °C during expiration.

Fig. 1 Temperature variation around the nostrils during inspiration and expiration. Thermograms of the nose during (a) inhalation, and (b) exhalation. The mean temperature of the ROM reaches 31.17 °C and 31.44 °C during inspiration and expiration, respectively. Color maps are depicted on the right side of the thermograms.

Download Full Size | PDF

To measure BR in the recorded thermal sequence, the nose [region of interest (ROI)] must be automatically identified in the first frame. Furthermore, a rough tracking in subsequent frames needs to be performed to compensate the motion of the subject or patient. In order to improve the signal-to-noise ratio (SNR), a second ROI, denominated region of measurement (ROM), was defined. Finally, the breathing waveforms can be extracted and used for computing the BR. Figure 2 illustrates the mentioned steps used to extract this vital parameter.

Fig. 2 Graphic representation of the main steps used to extract BR from thermal imaging. (1) Video sequence is recorded. (2) An automatic detection of the nose (ROI) is performed. Subsequently, a first rough tracking of this region is carried out. (3) Identification of ROM which encloses the nostrils. (4) The breathing waveform is extracted from the mean temperature of the ROM. (5) Breathing rate estimation.

Download Full Size | PDF

The steps of automatic detection of nose ROI, tracking of the ROI, ROM identification and BR extraction are described in detail in sections 2.1–2.4.

The presented approach was implemented in MATLAB (MATLAB 2014a, The MathWorks Inc., Natick, MA). Moreover, it was tested on a 64-Bit Windows 7 computer with a quad-core Intel® Core™ i5-3450 3.10 GHz processor, 16 GB RAM and a solid-state drive (SSD). The analysis of the data was done offline.

2.1. Automatic detection of ROI

Face segmentation was the first stage of the proposed approach. It was performed according to the following three steps.

Step 1: Face was segmented by using the multi-level Otsu’s method [26]. In 1979, Nobuyuki Otsu [27] introduced a clustering-based image thresholding approach, which permits to separate an image in two classes, background and foreground. This algorithm uses discriminant analysis to estimate an optimal threshold value (T^★) by minimizing the within-class variance ( $σ_{W}^{2}$ ) or analogously by maximizing the between-class variance ( $σ_{B}^{2}$ ). The latter is governed by the equation $T^{★} = \underset{1 \leq T < L}{argmax} {σ_{B}^{2} (T) = σ^{2} (T) - σ_{W}^{2} (T)},$ where σ² stands for the total variance, T is a threshold value and, lastly, L represents the gray levels. This equation can be further expanded to $σ_{B}^{2} (T) = ω_{1} (T) {[μ_{1} (T) - μ]}^{2} + ω_{2} (T) {[μ_{2} (T) - μ]}^{2} .$ Here, ω_i corresponds to the probabilities of the two classes (background and foreground), μ represents the mean intensity of the original image, and μ_i is the mean intensity of the respective class. As a matter of fact, 1 and 2 corresponds to classes C1 and C2, respectively.
This approach was extended by Liao and co-workers [26] and further optimized by other research groups [28] in order to allow multi-thresholding. Therefore, Eq. (1) and Eq. (2) are now given by
${T_{1}^{★}, T_{2}^{★}, \dots, T_{M - 1}^{★}} = \underset{1 \leq T_{1} < \dots < T_{M - 1} < L}{argmax} {σ_{B}^{2} (T_{1}, T_{2}, \dots, T_{M - 1})},$ and $σ_{B}^{2} (T) = \sum_{i = 1}^{M} ω_{k} (T) {[μ_{k} (T) - μ]}^{2} .$ Here, it was assumed that there are M − 1 thresholds (T₁, T₂,...,T_M−1) as well as optimal thresholds ( $T_{1}^{★}, T_{2}^{★}, \dots, T_{M - 1}^{★}$ ), which maximize ( $σ_{B}^{2}$ ). The M − 1 thresholds separate the image in M classes (C₁ for [ $1, \dots, T_{1}^{★}$ ], C₂ for [ $T_{1}^{★} + 1, \dots, T_{2}^{★}$ ], ..., C_M for [ $T_{M - 1}^{★} + 1, \dots, L$ ).
Step 2: After applying the multi-level Otsu’s method, remaining background noise was removed. Hence, it was considered that only the region with the largest area in the binary image coincides to the subject’s face.
Step 3: The next step aims at finding the chin contour. It was presented in [20].

In order to find the area around the nostrils in the segmented face [see Fig. 3(b)], facts from human anatomy and physiology were used. It is important to note that the thermograms of the segmented faces were converted to gray scale [ranges between 0.0 (black) and 1.0 (full intensity)] [Fig. 3(b)]. The idea behind this was to normalize the thermograms, which facilitates e.g. the selection of threshold values in later steps. The periorbital region (PR), more precisely the medial canthus region, illustrated in [Fig. 3(a)], is one of the warmest areas of the face as a result of its vascular anatomy [29]. Figure 3(c) illustrates an example. By modifying the scale’s range of Fig. 3(b) to [0.85, 1.0], it is possible to show this characteristic of the human physiology. Hence, to detect both PRs, we limited our search to window A [Fig. 3(a) and Fig. 3(d)]. Here, it was assumed that just the pixels with an intensity superior to a threshold of 0.95 correspond to the medial canthus regions [Fig. 3(d) and Fig. 3(e)]. An extended analysis of thermograms of differents subjects has demonstrated that just the pixels on these regions present an intensity higher than this threshold. Since the nostrils are below the PRs, the search window was limited again as demonstrated in Fig. 3(a) and Fig. 3(f). Finally, the canny edge detector was applied to Window B. Figure 3(g) shows the edges of the nose. In order to find the nostrils as well as to define the limits/size of ROI that will be posteriorly tracked, vertical P_V (x) and horizontal P_H(y) projections were calculated as follows,

P_{H} (y) = \sum_{x = 0}^{P - 1} I (x, y),

P_{V} (x) = \sum_{y = 0}^{Q - 1} I (x, y) .

Here, I(x,y) stands for the binary value at pixel (x,y) and, P and Q represent the width and height of the binary image, respectively. While Fig. 3(h) and Fig. 3(i) depict horizontal and vertical projections, Fig. 3(j) illustrates the limits of the detected ROI.

Fig. 3 Steps to find the ROI. (a) Illustration of proportions of the head, periorbital region - PR (blue dotted line window), and search windows A (black solid line) and B (green pattern made of crossing lines). (b) Thermogram of the segmented face. (c) Level and windowing adjustment of (b) to allow the visualization of periorbital regions. (d) Thermogram correspondent to search window A. (e) Dilated binary image containing the hot spots, which coincide with periorbital regions. (f) Thermogram correspondent to search window B. (g) Binary image containing the edges of (f) after applying the canny edge detector. (h) Horizontal projection P_H(y) of (g). The maximal value of P_H(y) is marked with an asterisk. (i) Vertical projection P_V (x) of (g). The two asterisks show the left and right extremities of the nose. (j) Dotted line rectangle enclosing the ROI that will be posteriorly tracked.

Download Full Size | PDF

2.2. Tracking of ROI

The tracking algorithm proposed by Mei et al. [30] was adopted to track the ROI detected in Section 2.1. The approach of this research is based on sparse representation, where the ℓ₁ -regularized least-squares problem is used to achieve sparsity. The sparse representation was integrated in a particle filter framework. The particle filter is composed of two models: a state transition model and an observation model, also denominated appearance model.

In the former, conditional density p(x_t|x_t−1) is used to estimate the dynamics of the target object between two successive frames. It points out the state transition probability. Here, x_t and x_t−1 represent the state variables, describing shape and location of the target, at time t and t − 1, respectively. Since the transition only depends on the previous state, the transition of the state variable x is considered to be a Markov chain.

The latter (observation model) can be considered as a target/background classifier. It estimates the similarity between target candidate and target model. The observation likelihood is defined by the conditional density p(z_t|x_t) at each time t. Here, x_t is the state variable and z_t denotes the observation. The maximum likelihood p(z_t|x_t) is given by

p (z_{t} | x_{t}) = \prod_{j = 1}^{n} 𝒩 ({\tilde{z}}_{t} (j); 0, σ^{2}),

where 𝒩 stands for the Gaussian distribution, σ² represents the variance, n is the number of pixels in the appearance model, j corresponds to the jth pixel, and lastly, z̃_t is the approximation residual of observation z_t by the target model.

As previously mentioned, this tracking algorithm uses the ℓ₁-regularized least-squares approach. It aims to minimize the absolute values of the residuals according to,

min_{c} {‖ B c - y ‖}_{2}^{2} + λ {‖ c ‖}_{1} .

Here, B consists of (1) a target template set and (2) trivial template sets; c is composed of both target and trivial coefficients; y is an 1D vector that corresponds to the target candidate; and λ is a regularization parameter. Additionally, |.|₁ and |.|₂ stand for the ℓ₁ and ℓ₂ norms, respectively.

2.3. Identification of ROM

In order to improve the SNR, a second and smaller ROI, denominated ROM, was defined (see Fig. 2). It was computed for each tracked ROI. The ROM was defined taking into consideration the edges of the nose [see Fig. 3(g)]. They were found by using the Canny edge detector. Forthwith, horizontal and vertical projections were calculated as defined in Eq. (5) and Eq. (6) [see Fig. 3(h) and Fig. 3(i)].

2.4. Extraction of breathing waveform and signal processing

For each frame the mean temperature value s̄(t) of the ROM was calculated according to

\bar{s} (t) = \frac{1}{m n} \sum_{i = 0}^{m - 1} \sum_{j = 0}^{n - 1} s (i, j, t),

where s(i, j, t) is the image temperature at pixel (i, j) and time point t, m describes the width of the ROI, an lastly, n represents its length. In summary, the breathing signal corresponds to the mean temperature of the ROM for all frames.

The signal was further preprocessed by applying a second order Butterworth band-pass filter with lower and upper 3 dB cutoff frequency of 0.1 Hz and 0.85 Hz, respectively.

To extract the instantaneous breathing frequencies, the algorithm introduced by Brüser et al. [31, 32] was implemented. These authors proposed a new algorithm that allows a robust estimation of local beat-to-beat intervals in physiological time series [33]. It consists in sliding an adaptive short analysis window w_i[v] across the signal s[n]. In this case, for each window location, the local breath-to-breath interval T_i is estimated. The length of the analysis window L + 1 was selected in order to contain at least two breaths, as given by

w_{i} [v] = s [n_{i} + v], v \in {- L / 2, \dots, L / 2},

with w_i[v] centered at n_i and L = 2T_maxf_s. In the latter, T_max denotes the expected upper breathing limit and f_s describes the sampling rate of the signal.

In order to estimate the local breath-to-breath interval T_i from each analysis window w_i[v], three estimators [(1) adaptive window autocorrelation - AC, (2) adaptive window average magnitude difference function - AAMDF and (3) maximum amplitude pairs - MAP] were independently computed. As a matter of fact, each estimator was calculated for the interval length m ∈ {m_min,...,m_max}, where m_min = T_minf_s and m_max = T_maxf_s.

Adaptive window autocorrelation - E_AC[m]
The adaptive window correlation introduced by Brüser and co-workers is estimated for all interval lengths (discrete lags) m according to,
$E_{AC} [m] = \frac{1}{m} \sum_{v = 0}^{m} w [v] \cdot w [v - m] .$ In brief, this estimator E_AC computes the correlation between m samples to the right w[v] and to the left w[v − m] of the analysis window center w[0].
Adaptive window average magnitude difference function - E_AMDF[m]
The average magnitude difference function estimator in line with autocorrelation was firstly applied to pitch tracking. It aims to find the absolute difference between samples. Brüser et al. [31, 32] adapted this method by using an adaptive window. This new approach is given by
$E_{AAMDF} [m] = {(\frac{1}{m} \sum_{v = 0}^{m} | w [v] - w [v - m] |)}^{- 1} .$
Maximum amplitude pairs - E_MAP[m]
The current estimator can be intended as an indirect peak detector, since it considers only the signal amplitude. E_MAP[m] is define as follows
$E_{MAP} [m] = max_{v \in {0, \dots, m}} (w [v] + w [v - m]) .$ For each lag m, the maximum amplitude of any two samples is computed. Therefore, E_MAP[m] reaches its maximum when two peaks, separated by m samples, are included in the analysis window.

To combine the three previous estimators (E_AC[m], E_AAMDF[m], E_MAP[m]), a Bayesian fusion method was applied. Hence, the conditional probability of m being the correct breath-to-breath interval given the three estimators P(m|E_AC, E_AAMDF, E_MAP) can be calculated according to

P (m | E_{AC}, E_{AAMDF}, E_{MAP}) \propto P (m | E_{AC}) \cdot P (m | E_{AAMDF}) \cdot P (m | E_{MAP}) .

As a matter of fact, the three estimators (E_AC[m], E_AAMDF[m], E_MAP[m]) can be interpreted as probability density functions (P(m|E_AC), P(m|E_AAMDF) and P(m|E_MAP)) as described in [33].

3. Experimental Setup

Thermography sequences were recorded by using a long wave infrared (LWIR) camera, VarioCAM® HD head 820S/30 mm (InfraTec GmbH, Dresden, Germany). It has an uncooled infrared microbolometer focal plane array (FPA) with a spatial resolution of 1024 × 768 pixels. This camera detects infrared wavelengths in the spectral range of 7.5 – 14 μm and presents a thermal sensitivity better than 0.05 K at 30 °C. In addition, infrared thermograms were acquired with a frame rate of 30 frames per second (fps).

Eleven healthy volunteers (4 females and 7 males), between the ages of 21 and 31 (25.36 years ± 2.94), participated in the current experiment. For the measurements, the camera was sat atop a tripod, placed approximately 2 m away from the subjects who sat on a chair [Fig. 4(a)]. The experiments were performed in a temperature-controlled environment (room temperature and humidity were approximately 22 °C and 50%, respectively).

Fig. 4 (a) Experimental Setup (1 - body attached data recording system SOMNOlab 2 (Weinmann GmbH, Hamburg, Germany), 2 - Thermographic camera VarioCAM® HD head 820S/30 mm (InfraTec GmbH, Dresden, Germany)). (b) Small sequence of the breathing waveform displayed in the monitor to simplify the subjects’ task.

Download Full Size | PDF

The study protocol consisted of three phases A to C. In phase A, the subjects were instructed to stay still and to breathe normally. A three minutes recording (Sequence 1) was performed. During phase B, the subjects were advised to breathe normally and to perform some movements during the three minutes acquisition (Sequence 2) in order to test the robustness of our approach against motion. In the last phase C, they were asked to follow a sequence containing different breathing patterns (eupnea, tachypnea, apnea, deep breathing (Kussmaul breathing) and Cheyne-Stokes respiration). It is depicted in Fig. 5. In order to simplify this task, breathing patterns were shown on a monitor. The subjects followed the displayed breathing waveform as demonstrated in the example of Fig. 4(b). The sequence recorded corresponds to sequence 3.

Fig. 5 Intended breathing pattern sequence used during phase C.

Download Full Size | PDF

Thoracic effort (piezoplethysmography) was the ground truth chosen to validate our results. To measure this parameter the data recording system SOMNOlab 2 (Weinmann GmbH, Hamburg, Germany) was utilized.

4. Results

During phase A, the subjects breathed normally and stayed still according to the study protocol. Table 1 shows the performance of our algorithm for each subject. On average, 97.52% of each breathing waveform did not contain any artifact and was usable for BR estimation (in 10 out of 11 subjects coverage was 100%). A mean correlation (Corr.) between BR estimated with IRT and BR obtained with piezoplethysmography reached 0.968 ± 0.030. Moreover, the mean absolute BR error (Ē) was 0.33 bpm and the spread of the error, calculated using the 90^th percentile of the errors (E₉₀), averaged 0.71 bpm. Figure 6(a) depicts a Bland-Altman plot comparing both measurement techniques, IRT and ground truth, for subject S9. The estimated mean difference is 0.004 bpm and the limits of agreement range from −1.2 bpm to 1.2 bpm.

Fig. 6 Bland-Altman plots for subjects S9 (phase A), S3 (phase B) and S10 (phase C), respectively. They compare BR obtained with IRT (BR_IRT) with reference BR (BR_Ref). (a) This first plot, correspondent to phase A, shows a bias of 0.004 bpm and limits of agreement of −1.2 bpm to 1.2 bpm. (b) In the second example (subject 3, phase B) the bias is −0.047 and the limits of agreement range from −1.4 bpm to 1.3 bpm. (c) The last Bland-Altman plot (subject 10, phase C) presents a mean difference of 0.12. The 95% limits of agreement vary between −3.7 and 3.9 bpm. As a matter of fact, solid line indicates the bias/mean difference and 95% of limits of agreement are plotted in dashed lines.

Download Full Size | PDF

Table 1. Performance of the proposed algorithm for phase A

View Table | View all tables in this article

During phase B, subjects performed some head movements while breathing normally. Our approach was again able to estimate the respiratory rate of all subjects. In this case 99.36% of the data was free from artifacts, and thus, suitable for BR extraction. In Table 2 the results of our algorithm are presented. It shows a mean correlation of 0.940 ± 0.032. Furthermore, an average breathing rate error of 0.55 bpm was achieved. The spread of the error reaches, in turn, 1.23 bpm. The Bland-Altman plot for subject S3 is represented in Fig. 6(b). In this example, the bias is −0.047 and the limits of agreement vary between −1.4 bpm and 1.3 bpm. Figure 7(a) represents a thermogram of sequence 2 for the same volunteer as well as the respective breathing waveform. It is contained in Visualization 1, which aims to illustrate the capacity of the algorithm to follow a moving ROI and still be able to extract the breathing waveform.

Fig. 7 (a) Thermogram of sequence 2 for subject S3. A short animation illustrating the performance of the tracking algorithm is shown in Visualization 1. (b) Thermogram of sequence 3 for subject S10. A short animation illustrating the performance of the algorithm is shown in Visualization 2.

Download Full Size | PDF

Table 2. Performance of the proposed algorithm for phase B

View Table | View all tables in this article

During phase C, subjects simulated different breathing patterns (see Fig. 5). Table 3 demonstrates that all data was free from artifacts, and therefore, usable for breathing rate estimation. The mean correlation corresponds to 0.974 ± 0.017. Lastly, the BR error and the spread of the error averaged 0.96 bpm and 1.51 bpm, respectively. A Bland-Altman plot for subject S10 is depicted in Fig. 6(c). Here, the mean difference averages 0.12 and the 95% limits of agreement range from −3.7 to 3.9 bpm. Figure 7(b) depicts a thermogram of sequence 3 for the same candidate as well as the breathing waveform. It is a single frame of the animation that can be visualized in Visualization 2. It illustrates the capacity of the algorithm to be sensitive to different breathing patterns. Figure 8 shows the spectrogram of the third sequence for subject S4. In addition, it depicts the BR estimated with the proposed approach (solid line) as well as the BR obtained with the ground truth (dashed line). Three breathing waveforms segments (A, B and C), corresponding to eupnea, apnea and tachypnea, are also represented. Again, solid line and dashed line represent IRT signal and ground truth signal, respectively.

Fig. 8 Top: Spectrogram of the third IRT sequence of subject S4. Middle: Estimated BR (solid line - IRT, dashed line - piezoplethysmography). Bottom: Three 30-seconds segments of IRT (solid line) and piezoplethysmography waveforms (dashed line). The segments A, B and C correspond to eupnea, apnea, tachypnea, respectively.

Download Full Size | PDF

Table 3. Performance of the proposed algorithm for phase C

View Table | View all tables in this article

5. Discussion

Breathing rate is one of the primary vital signs. It can be an important indicator of treatment effects, procedural complications or clinical deterioration. Monitoring of respiratory function is fundamental. However, studies have demonstrated that it is one of the most frequently undocumented parameters [9]. This intensifies the need for reliable unobtrusive monitoring alternatives.

The present scientific paper introduces a new approach that allows to remotely monitor breathing rate (BR) by using thermography. To evaluate the performance of the developed algorithm, and thus, the feasibility of IRT as a valid noncontact measurement technique, a study in eleven healthy participants was conducted. It focused on analyzing the robustness of our method under challenging conditions: (a) motion (b) variable breathing rate and (c) respiratory disorders. Therefore, three sequences for each subject were acquired.

During the first phase the aim was to investigate and compare the performance of the algorithm without any movement artifacts. Table 1 shows an excellent agreement between both methods (infrared thermography and piezoplethysmography) with a correlation coefficient of 0.97 ± 0.03 and a mean absolute BR error that did not exceed 0.33 bpm. Lewis and coworkers [22] published in 2011 a paper where they used infrared video thermography to estimate breathing rate and relative tidal volume. Using a SC-6000 [320 × 240 pixel resolution and sensitivity of 0.08 °C (Indigo System Inc., Goleta, CA)] and a TVS-700 [640 × 510 pixel resolution and sensitivity of 0.02 °C (FLIR Inc., Santa Barbara, CA)] infrared imaging camera they obtained similar correlations, 0.95 ± 0.05 (N=6) and 0.98 ± 0.02 (N = 12), respectively. Three 2 minutes recordings were acquired, in which the subjects were advised to breathe through the nose (a) normally/spontaneously, (b) slowly and deeply, and (c) rapidly. In addition, they were instructed to minimize head movements. In 2010 Fei et al. [12] published a paper where they also used IRT [IRT camera FLIR SC6000, 640 × 542 pixel resolution and sensitivity of 0.025 °C (FLIR Inc., Santa Barbara, CA)] to recover the breathing signal from the subjects’ nostrils. They recorded approximately 3 minutes thermal clips for twenty candidates. To validate their method, a measure denominated CAND (Complement of the Absolute Normalized Difference) was proposed. It is given by

CAND = 1 = \frac{| \bar{{BR}_{T}} - \bar{{BR}_{G}} |}{\bar{{BR}_{G}}},

where

\bar{{BR}_{T}}

is the mean BR obtained with IRT and

\bar{{BR}_{G}}

corresponds to the mean BR obtained with the ground truth. The mean CAND for the experiment reached approximately 98%. In summary, the three approaches present outstanding results as demonstrated by the high correlations between IRT and gold standard. Figure 6(a) shows a Bland-Altman plot of the relation between IRT and piezoplethysmography correspondent to volunteer S9. The estimated bias (0.004 bpm) indicates that in this case the two methods have very similar results on average. In addition, a good agreement between both methods can be observed here.

During phase B, subjects carried out arbitrary head movements. Since we always wanted to have the nose in the camera’s field of view (FOV) the subjects tried to avoid strong rotations of the head (down-, left- and rightwards). Visualization 1 demonstrates that despite rapid movements the tracking algorithm was able to efficiently follow the ROI. Table 2 shows again an excellent concordance between thermal imaging and piezoplethysmography. The correlation coefficient was 0.94 ± 0.03 with a mean absolute BR error lower than 0.55 bpm. In Fig. 6(b) a Bland-Altman plot comparing both methods for candidate S3 is presented. Also here IRT and ground truth have similar results on average and an excellent agreement. To the best of our knowledge, there is no further published study that examines the effect of movements on the algorithm’s performance. Indeed, Lewis et al. [22] considered that as a limitation of their and other works, since they applied tracking algorithms but did not validate them.

In order to examine (simulated) clinically relevant scenarios within our study, a broad range of breathing patterns, such as eupnea, tachypnea, apnea, Kussmaul breathing and Cheyne-Stokes respiration, were simulated. Also here we obtained promising results (see Visualization 2). The correlation between IRT and gold standard is higher than 0.97 and the mean absolute breathing rate error was 0.96 bpm (E₉₀=1.51 bpm). A Bland-Altman plot correlating both methods is depicted in Fig. 6(c). It corresponds to volunteer S10. The obtained mean difference (0.12 bpm) evidences that the two methods have similar results on average. Furthermore, a good agreement between IRT and gold standard can be here observed. Figure 8 illustrates, in turn, the spectrogram of the IRT waveform for candidate S4 as an example. Frequencies with higher amplitude, which coincide with respiratory rate, are depicted in dark red. This visual representation describes the variability of BR caused by the different breathing patterns very well (compare with Fig. 5). The second plot shows the BR estimated from both methods. An excellent concordance between the signals is noticeable. Additionally, Fig. 8 also shows three 30-seconds segments (A - eupnea, B - apnea, C - tachypnea) of sequence 3. These three intervals are examples of the excellent agreement between signals of the two modalities, IRT and piezoplethysmography.

Figure 8 (middle) together with Table 3 show a high mean absolute BR error when compared with Table 1 and 2. These errors are also visible in Fig. 6(c), which correspond to the outliers at approximately 6 bpm and 20 bpm. In phase C abrupt changes in breathing rate were simulated, therefore, small delays between both signals in these transitions caused errors of about 15 bpm, which negatively influenced the final mean absolute error given in Table 3. The disparity between signals occurred in the transitions breathing/apnea, apnea/breathing, eupnea/tachypnea and tachypnea/eupnea as you can confirm in the example given in Fig. 8 (middle). It further indicates that the error between both methods outside the transitions is quite small. Furthermore, we are comparing methods that measure the same vital parameter but are based on different measurement principles. This might also contribute to an increase of the errors between both methods. However, we strongly believe that our approach might be appropriated for clinical purposes. It is able to correctly extract breathing rate in different scenarios (e.g. tachypnea, apnea, Kussmaul breathing, and Cheyne-stokes respiration). Therefore, the algorithm used to estimate BR may be a good alternative for nonstationary signals, including breathing rate and heart rate. Moreover, since the delays between both methods in the transitions periods are quite small, we did not think that it can be a great issue in clinical monitoring.

In contrast to phase C, in phase A and B a 100% coverage was not achieved. The infrared camera used in this study presents a Non Uniformity Correction (NUC) system. NUC permits to homogenize the image, i.e. it corrects the low temperature deviations of the detector elements. During this correction, performed usually within a fixed time interval, the live display and image acquisition is briefly interrupted. Normally, NUC has no impact on our algorithm. However, during the video acquisition of two subjects (S2 phase A, S1 phase B) the camera made within a significant period of time successive corrections. Since too much information was lost, the signal presents a lot of artifacts. Therefore, we did not consider these intervals for data evaluation.

There is a wide range of research papers in the literature that focus on other technologies (e.g. Doppler radar [1] and visible/near infrared image sensors [11]) to measure breathing activity. Together with IRT, these two technologies have demonstrated a great potential to detect this important vital parameter. However, they present some disadvantages in comparison with thermal imaging. For example, Doppler radar is extremely sensitive to motion artifacts and to the presence of multiple subjects. Visible/near infrared imaging systems need, in turn, an illumination source.

6. Conclusions

In the last decades there is a need for unobtrusive monitoring of vital signs. A wide range of research groups have focused their work on camera-based monitoring alternatives. Visible, near infrared and long-wave infrared are some examples of the technologies that are investigated.

This scientific paper introduces a new robust and feasible approach which uses infrared imaging for monitoring of breathing function. Our method includes a tracking algorithm to compensate the effect of head movements on the IRT waveform extraction. Moreover, it uses a powerful algorithm for BR estimation which has demonstrated to be suitable for nonstationary signals. The results show that infrared thermography is able to accurately estimate breathing rate under challenging conditions such as motion and possible respiratory disorders.

In a further step, reliability and feasibility of this approach in a clinical scenario should be taken into consideration. Equally important is the impact of room temperature on the temperature measured around the nostrils. It is an interesting issue that must be investigated in future studies. In addition, to obtain a more robust extraction of this vital parameter, we are currently developing a method that considers temperature variation around the nostrils and mouth as well as the movement of the shoulders during inhalation and exhalation. Apart from nasal breathing, open-mouth respiration is also normal in adult humans. Shoulders and thorax movement are other characteristics that are prominent during this process. They can be useful if neither nose nor mouth are in the field of view of the infrared camera. Therefore, it would be propitious to fuse all sources to permit a continuous robust and feasible monitoring of this crucial vital parameter.

Acknowledgments

C. B. Pereira wishes to acknowledge FCT (Foundation for Science and Technology in Portugal) for her PhD grant SFRH / BD / 84357 / 2012.

References and links

1. A. D. Droitcour, T. B. Seto, B.-K. Park, S. Yamada, A. Vergara, C. El Hourani, T. Shing, A. Yuen, V. M. Lubecke, and O. Boric-Lubecke, “Non-contact respiratory rate measurement validation for hospitalized patients,” in Proceedings of IEEE Engineering in Medicine and Biology Society (IEEE, 2009), pp. 4812–4815.

2. M. A. Cretikos, R. Bellomo, K. Hillman, J. Chen, S. Finfer, and A. Flabouris, “Respiratory rate: the neglected vital sign,” Med. J. Aust. 188(11), 657–659 (2008). [PubMed]

3. J. Fei and I. Pavlidis, “Virtual Thermistor,” in Proceedings of IEEE Engineering in Medicine and Biology Society (IEEE, 2007), pp. 250–253.

4. R. K. Modak, Anesthesiology Keywords Review (Lippincott Williams & Wilkins, 2013).

5. R. Riha, “Diagnostic approaches to respiratory sleep disorders,” J. Thorac. Dis. 7(8), 1373–1384 (2015).

6. T. Duning, M. Deppe, E. Brand, J. Stypmann, C. Becht, A. Heidbreder, and P. Young, “Brainstem involvement as a cause of central sleep apnea: pattern of microstructural cerebral damage in patients with cerebral microangiopathy,” PLoS One 8(4), e60304 (2013). [CrossRef] [PubMed]

7. M. Siccoli, P. Valko, D. Hermann, and C. Bassetti, “Central periodic breathing during sleep in 74 patients with acute ischemic stroke - neurogenic and cardiogenic factors,” J. Neurol. 255(11), 1687–1692 (2008). [CrossRef] [PubMed]

8. M. Cretikos, J. Chen, K. Hillman, R. Bellomo, S. Finfer, and A. Flabouris, “The objective medical emergency team activation criteria: A case-control study,” Resuscitation 73(1), 62–72 (2007). [CrossRef] [PubMed]

9. F. McGain, M. Cretikos, D. Jones, S. V. Dyk, M. D. Buist, H. Opdam, V. Pellegrino, M. S. Robertson, and R. Bellomo, “Documentation of clinical review and vital signs after major surgery,” Med. J. Aust. 189(7), 380–383 (2008). [PubMed]

10. A. K. Abbas, K. Heimann, K. Jergus, T. Orlikowsky, and S. Leonhardt, “Neonatal non-contact respiratory monitoring based on real-time infrared thermography,” Biomed. Eng. Online 10(93), 1–17 (2011). [CrossRef]

11. M. Bartula, T. Tigges, and J. Muehlsteff, “Camera-based system for contactless monitoring of respiration,” in Proceedings of IEEE Engineering in Medicine and Biology Society (IEEE, 2013), pp. 2672–2675.

12. J. Fei and I. Pavlidis, “Thermistor at a distance: unobtrusive measurement of breathing,” IEEE Trans. Biomed. Eng. 57(4), 988–998 (2010). [CrossRef]

13. O. Mimoz, T. Benard, A. Gaucher, D. Frasca, and B. Debaene, “Accuracy of respiratory rate monitoring using a non-invasive acoustic method after general anaesthesia,” Br. J. Anaesth. 108(5), 872–875 (2012). [CrossRef] [PubMed]

14. J. H. Atkins and J. E. Mandel, “Performance of Masimo rainbow acoustic monitoring for tracking changing respiratory rates under laryngeal mask airway general anesthesia for surgical procedures in the operating room: a prospective observational study,” Anesth. Analg. 119(6), 1307–1314 (2014). [CrossRef] [PubMed]

15. Y. Guechi, A. Pichot, D. Frasca, F. Rayeh-Pelardy, J.-Y. Lardeur, and O. Mimoz, “Assessment of noninvasive acoustic respiration rate monitoring in patients admitted to an Emergency Department for drug or alcoholic poisoning,” J. Clin. Monit. Comput. (2015) (Epub ahead of print). [CrossRef] [PubMed]

16. S. Leonhardt and B. Lachmann, “Electrical impedance tomography: the holy grail of ventilation and perfusion monitoring?” Intensive Care Med. 38(12), 1917–1929 (2012). [CrossRef] [PubMed]

17. M. Czaplik, I. Biener, R. Dembinski, P. Pelosi, T. Soodt, W. Schroeder, S. Leonhardt, G. Marx, R. Rossaint, and J. Bickenbach, “Analysis of regional compliance in a porcine model of acute lung injury,” Respir. Physiol. Neurobiol. 184(1), 16–26 (2012). [CrossRef] [PubMed]

18. R. Murthy, I. Pavlidis, and P. Tsiamyrtzis, “Touchless monitoring of breathing function,” in Proceedings of IEEE Engineering in Medicine and Biology Society (IEEE, 2004), pp. 1196–1199.

19. R. G. Soto, E. S. Fu, H. Vila, and R. V. Miguel, “Capnography accurately detects apnea during monitored anesthesia care,” Anesth. Analg. 99(2), 379–382 (2004). [CrossRef] [PubMed]

20. C. B. Pereira, X. Yu, V. Blazek, and S. Leonhardt, “Robust remote monitoring of breathing function by using infrared thermography,” in Proceedings of IEEE Engineering in Medicine and Biology Society (IEEE, 2015), pp. 4250–4253.

21. M. Garbey, N. Sun, A. Merla, and I. Pavlidis, “Contact-free measurement of cardiac pulse based on the analysis of thermal imagery,” IEEE Trans. Biomed. Eng. 54(8), 1418–1426 (2007). [CrossRef] [PubMed]

22. G. F. Lewis, R. G. Gatto, and S. W. Porges, “A novel method for extracting respiration rate and relative tidal volume from infrared thermography,” Psychophysiology 48(7), 877–887 (2011). [CrossRef] [PubMed]

23. K. Heimann, K. Jergus, A. K. Abbas, N. Heussen, S. Leonhardt, and T. Orlikowsky, “Infrared thermography for detailed registration of thermoregulation in premature infants,” J. Perinat. Med. 41(5), 613–620 (2013). [PubMed]

24. C. B. Pereira, M. Czaplik, N. Blanik, R. Rossaint, V. Blazek, and S. Leonhardt, “Contact-free monitoring of circulation and perfusion dynamics based on the analysis of thermal imagery,” Biomed. Opt. Express 5(4), 1075–1089 (2014). [CrossRef] [PubMed]

25. R. Chauvin, M. Hamel, S. Brière, F. Ferland, F. Grondin, D. Létourneau, M. Tousignant, and F. Michaud, “Contact-free respiration rate monitoring using a pan-tilt thermal camera for stationary bike telerehabilitation sessions,” IEEE Syst. J. PP(99), 1–10 (2014). [CrossRef]

26. P. S. Liao, T. S. Chen, and P. C. Chung, “A fast algorithm for multi-level thresholding,” J. Inf. Sci. Eng. 17(5), 713–727 (2001).

27. N. Otsu, “A Threshold Selection Method from Gray-Level Histograms,” IEEE Trans. Syst. Man Cybern. 9(1), 62–66 (1979). [CrossRef]

28. D.-Y. Huang and C.-H. Wang, “Optimal multi-level thresholding using a two-stage Otsu optimization approach,” Pattern Recognit. Lett. 30(3), 275–284 (2009). [CrossRef]

29. Z. Zhu, J. Fei, and I. Pavlidis, “Tracking human breath in infrared imaging,” in Proceedings of IEEE Symposium on Bioinformatics and Bioengineering (IEEE, 2005), pp.227–231.

30. X. Mei and H. Ling, “Robust Visual Tracking and Vehicle Classification via Sparse Representation,” IEEE Trans. Pattern Anal. Mach. Intell. 33(11), 2259–2272 (2011). [CrossRef] [PubMed]

31. C. Brüser, S. Winter, and S. Leonhardt, “How speech processing can help with beat-to-beat heart rate estimation in ballistocardiograms,” in Proceedings of IEEE Medical Measurements and Applications (IEEE, 2013), pp.12–16.

32. C. Brüser, S. Winter, and S. Leonhardt, “Robust inter-beat interval estimation in cardiac vibration signals,” Phys. Meas. 34(2), 123–138 (2013). [CrossRef]

33. C. Brüser, Sensing and Signal Processing for Unobtrusive Cardiac Monitoring Utilizing Ballistocardiography (Shaker Verlag, 2014).

Name	Description
Visualization 1: MP4 (6834 KB)	Animation illustrating the performance of the tracking algorithm.
Visualization 2: MP4 (9919 KB)	Animation illustrating the performance of the algorithm during dynamic breathing rate.

Subject	Gender	Duration (min)	Coverage (%)	Correlation		BR_abs	Error (bpm)

				Corr.	^*CI₉₅	Ē	E₉₀
S1	F	3	100	0.979	[0.975; 0.982]	0.32	0.65
S2	M	2.18	72.67	0.901	[0.883; 0.916]	0.56	1.13
S3	M	3	100	0.986	[0.985; 0.987]	0.21	0.53
S4	M	3	100	0.986	[0.985; 0.987]	0.26	0.56
S5	M	3	100	0.990	[0.989; 0.991]	0.18	0.42
S6	F	3	100	0.924	[0.918; 0.930]	0.64	1.18
S7	M	3	100	0.978	[0.976; 0.980]	0.24	0.54
S8	F	3	100	0.976	[0.974; 0.978]	0.35	0.75
S9	M	3	100	0.952	[0.948; 0.956]	0.44	0.98
S10	M	3	100	0.990	[0.989; 0.991]	0.22	0.50
S11	F	3	100	0.986	[0.985; 0.988]	0.25	0.53

Mean		2.93	97.52	0.968		0.33	0.71

Subject	Gender	Duration (min)	Coverage (%)	Correlation		BR_abs Error (bpm)

				Corr.	^*CI₉₅	Ē	E₉₀
S1	F	2.79	93	0.892	[0.881; 0.903]	0.72	1.86
S2	M	3	100	0.941	[0.934; 0.947]	0.62	1.34
S3	M	3	100	0.961	[0.958; 0.964]	0.45	0.78
S4	M	3	100	0.962	[0.959; 0.965]	0.43	0.89
S5	M	3	100	0.912	[0.903; 0.920]	0.80	1.55
S6	F	3	100	0.880	[0.868; 0.890]	0.76	1.85
S7	M	3	100	0.980	[0.979; 0.982]	0.31	0.75
S8	F	3	100	0.957	[0.953; 0.960]	0.50	1.21
S9	M	3	100	0.965	[0.962; 0.967]	0.43	1.05
S10	M	3	100	0.955	[0.952; 0.959]	0.47	1.18
S11	F	3	100	0.931	[0.925; 0.936]	0.61	1.11

Mean		2.98	99.36	0.940		0.55	1.23

Subject	Gender	Duration (min)	Coverage (%)	Correlation		BR_abs Error (bpm)

				Corr.	^*CI₉₅	Ē	E₉₀
S1	F	10.43	100	0.963	[0.957; 0.968]	1.36	2.17
S2	M	10.03	100	0.938	[0.932; 0.944]	1.45	1.96
S3	M	10.26	100	0.982	[0.982; 0.983]	0.72	1.27
S4	M	10.15	100	0.949	[0.942; 0.955]	1.37	1.64
S5	M	10.24	100	0.983	[0.982; 0.984]	0.64	0.71
S6	F	9.97	100	0.988	[0.987; 0.989]	0.82	1.65
S7	M	10.21	100	0.992	[0.992; 0.993]	0.58	1.05
S8	F	9.71	100	0.978	[0.976; 0.979]	0.98	1.58
S9	M	10.38	100	0.978	[0.976; 0.979]	1.12	2.05
S10	M	9.94	100	0.990	[0.990; 0.991]	0.65	1.07
S11	F	10.30	100	0.976	[0.974; 0.977]	0.87	1.41

Mean		10.15	100	0.974		0.96	1.51

Subject	Gender	Duration (min)	Coverage (%)	Correlation		BR_abs	Error (bpm)

				Corr.	^*CI₉₅	Ē	E₉₀
S1	F	3	100	0.979	[0.975; 0.982]	0.32	0.65
S2	M	2.18	72.67	0.901	[0.883; 0.916]	0.56	1.13
S3	M	3	100	0.986	[0.985; 0.987]	0.21	0.53
S4	M	3	100	0.986	[0.985; 0.987]	0.26	0.56
S5	M	3	100	0.990	[0.989; 0.991]	0.18	0.42
S6	F	3	100	0.924	[0.918; 0.930]	0.64	1.18
S7	M	3	100	0.978	[0.976; 0.980]	0.24	0.54
S8	F	3	100	0.976	[0.974; 0.978]	0.35	0.75
S9	M	3	100	0.952	[0.948; 0.956]	0.44	0.98
S10	M	3	100	0.990	[0.989; 0.991]	0.22	0.50
S11	F	3	100	0.986	[0.985; 0.988]	0.25	0.53

Mean		2.93	97.52	0.968		0.33	0.71

Remote monitoring of breathing dynamics using infrared thermography

Abstract

1. Introduction

2. Methodology

2.1. Automatic detection of ROI

2.2. Tracking of ROI

2.3. Identification of ROM

2.4. Extraction of breathing waveform and signal processing

3. Experimental Setup

4. Results

5. Discussion

6. Conclusions

Acknowledgments

References and links

Supplementary Material (2)

Cited By

Figures (8)

Tables (3)

Equations (15)

Biomedical Optics Express