figure table

# Development of a blood oxygenation phantom for photoacoustic tomography
combined with online pO${}_{\textbf{2}}$ detection and flow spectrometry

Marcel Gehrung  Cancer Research UK Cambridge Institute, Li Ka-Shing Centre,
Cambridge, UK  Department of Physics, University of Cambridge, Cambridge, UK
Sarah E. Bohndiek  Cancer Research UK Cambridge Institute, Li Ka-Shing Centre,
Cambridge, UK  Department of Physics, University of Cambridge, Cambridge, UK
Joanna Brunker  Cancer Research UK Cambridge Institute, Li Ka-Shing Centre,
Cambridge, UK  Department of Physics, University of Cambridge, Cambridge, UK

###### Abstract

Photoacoustic tomography (PAT) is intrinsically sensitive to blood oxygen
saturation (sO2) in vivo. However, making accurate sO2 measurements without
knowledge of tissue- and instrumentation-related correction factors is
extremely challenging. We have developed a low-cost flow phantom system to
facilitate validation of photoacoustic tomography systems. The phantom is
composed of a flow circuit, which is partially embedded within a tissue
mimicking phantom, with independent sensors providing online monitoring of the
optical absorption spectrum and partial pressure of oxygen in the tube. We
first establish the flow phantom using two small molecule dyes that are
frequently used for photoacoustic imaging: methylene blue (MB) and indocyanine
green (ICG). We then demonstrate the potential of the phantom for evaluating
sO2 using chemical oxygenation and deoxygenation of blood in the phantom.
Using this dynamic assessment of the photoacoustic sO2 measurement in phantoms
in relation to a ground truth, we explore the influence of multispectral
processing and spectral coloring on accurate assessment of sO2. Future studies
could exploit this low-cost dynamic flow phantom to validate fluence
correction algorithms and explore additional blood parameters such as pH, and
also absorptive and other properties of different fluids.

###### keywords:

blood oxygenation, flow, phantom, photoacoustic tomography

*Joanna Brunker<EMAIL_ADDRESS>
##  1 Introduction

Photoacoustic tomography (PAT) exploits optically generated ultrasound to
provide images that combine the high contrast and spectral specificity of
optical imaging with the high spatial resolution of ultrasound. In particular,
PAT has been widely used to image blood hemoglobin concentration and
oxygenation, which have the potential to inform on a range of
pathophysiologies, from tumour aggressiveness [1] and treatment response [2,
3] to intestinal inflammation associated with Crohn’s disease [4] and colitis
[5]. The derivation of such images is usually based on spectral unmixing to
resolve the differential absorption contributions of oxy- (HbO2) and deoxy-
hemoglobin (Hb). Total hemoglobin concentration (THb) is typically taken as
the sum of the contributions to the photoacoustic signal
$P(\mathbf{r},\lambda)$ from HbO2 and Hb, while hemoglobin oxygenation (sO2)
is taken as the ratio of HbO2 to THb.

Unfortunately, the assessment of HbO2 and Hb content from photoacoustic data
is not trivial. Estimation of these chromophore concentrations from images
taken at multiple wavelengths is commonly achieved using least-squares fitting
of reference HbO2 and Hb spectra [6], which are tabulated in the literature
[7, 8, 9]; however, these spectra are recorded under in vitro conditions thus
can vary substantially from those experienced within an in vivo study.
Moreover, the measured photoacoustic signal $P(\mathbf{r},\lambda)$ is not
directly proportional to the absorbed energy density $H(\mathbf{r},\lambda)$.
$H(\mathbf{r},\lambda)$ is the product of the light fluence $\Phi$, which
itself varies as a function of $\mathbf{r}$ and $\lambda$, as well as
$\mu_{a}(\mathbf{r},\lambda)$ and the reduced scattering coefficient
$\mu^{\prime}_{s}(\mathbf{r},\lambda)$:

$H(\mathbf{r},\lambda)=\mu_{a}(\mathbf{r},\lambda)\Phi[\mathbf{r},\lambda,\mu_{a}(\mathbf{r},\lambda),\mu^{\prime}_{s}(\mathbf{r},\lambda)].$
(1)

This co-dependence of $H(\mathbf{r},\lambda)$ on both absorption and light
fluence leads to an effect known as “spectral coloring”, where variations in
local fluence bias the measured optical absorption distribution [10, 11].
Despite this challenge, the amplitude and temporal profiles of photoacoustic
waves have been shown to be linearly dependent on the blood sO2 in vitro [12].
In theory, only two wavelengths are necessary to calculate two concentrations,
but practical limitations mean the choice and number of wavelengths becomes
critical for accurate sO2 measurements [13]. Alternative approaches beyond the
standard linear unmixing model have been shown to improve sO2 measurement
accuracy, in particular by accounting for the absorbed energy distribution
through: internal irradiation [14]; diffusion theory modelling [15, 16, 17,
18, 19]; Monte Carlo simulations [20]; model-based iterative minimisation [21,
22, 23]; and linear superposition of reference fluence base spectra [24].

Several studies have explored ways to validate photoacoustic images using
phantoms with well-characterised optical and acoustic properties [25, 26].
Validation of photoacoustically measured sO2 is possible using CO-oximetry
[12, 15, 27], pulse oximetry [28] and blood-gas analysis [29]; correlation
between pO2 and sO2 has also been used to study oxygen-hemoglobin binding
[30]. However, most studies have worked with static blood samples with limited
control on the sO2. Allowing blood to flow in a circuit [27] and with the
ability to vary the sO2 provides a more versatile platform for investigating
PAT oxygenation measurements.

In this work, we created a low-cost flow phantom system with online monitoring
to facilitate validation of PAT systems. The phantom was tested by circulating
different concentrations of methylene blue (MB) and indocyanine green (ICG).
The phantom was then applied to explore the accuracy of sO2 assessment using
PAT, with sO2 values ranging from $0\text{\,}\mathrm{\char 37\relax}$ to
$100\text{\,}\mathrm{\char 37\relax}$. We studied the impact of evaluating sO2
using different reference spectra for HbO2 and Hb, and also the effect of
spectral coloring, showing discrepancies of up to $60\text{\,}\mathrm{\char
37\relax}$ between the actual and measured sO2. Our findings highlight the
importance of careful choice of spectra for unmixing and the development of
fluence correction models to improve the biological relevance of sO2
measurements derived from PAT images.

##  2 Methods

###  2.1 Flow system

The flow system (Fig. 1) enables fluids to be circulated within a vessel-
mimicking tube embedded in a tissue-mimicking agar phantom placed in the
chamber of a photoacoustic imaging system. Online (and offline) flow
spectrometry and pO2 detection provide independent validation of the spectral
measurements made by the imaging system. Details of the tissue-mimicking
phantom, the photoacoustic imaging system, the spectrometric validation, the
pO2 detection, and the spectral analysis are described in the following
subsections.

Figure 1: Overview of the flow system. (1) Injection site for introducing
oxygenated blood (or other fluids) into the flow system, and for subsequently
deoxygenating the blood using sodium hydrosulfite delivered via the syringe
driver (MKCB2159V, Harvard); (2) spectra are recorded using a light source
(Avalight-HAL-S-Mini, Avantes) and spectrometer (AvaSpec-ULS2048-USB2-VA-50,
Avantes) as the blood passes through a flow cell (170700-0.5-40, Hellma
Analytics); (3) needle probes (NX-BF/O/E, Oxford Optronix) measure the
temperature and partial pressure of oxygen (pO2) before (3a) and after (3b)
the blood passes through an agar phantom immersed in the photoacoustic imaging
system (MSOT inVision 256-TF, iThera Medical) (4); (5) a touch-screen monitor
(OxyLite Pro, Oxford Optronix) displays temperature and oxygen data; (6) these
data are downloaded via an Arduino UNO and read in MATLAB on a laptop, which
also records the spectrometer readings via AvaSoft software; (7) a peristaltic
pump (CTP100, Fisher Scientific) provides blood circulation.

####  2.1.1 Tissue-mimicking phantom

The tissue-mimicking phantom used for photoacoustic imaging comprised either
static or flowing fluids within a tube (“vessel”) embedded within an agar
(“tissue”) cylinder, $19\text{\,}\mathrm{mm}$ in diameter. Agar phantoms were
prepared by heating a solution of 1.5% w/v agar (05039, Fluka) in water, and
then adding 2.1% v/v Intralipid (I141, Sigma-Aldrich), before pouring into an
open $20\text{\,}\mathrm{mL}$ syringe (with the injection end removed), in the
centre of which was positioned a needle supporting the tubing vertically. The
solidified agar phantom had a scattering coefficient of
$5\text{\,}{\mathrm{cm}}^{-1}$ due to the Intralipid, but the absorption
coefficient was assumed to be negligible.

To provide optical absorption for assessment of spectral coloring, nigrosin
(198285, Sigma-Aldrich) was also added to the tissue-mimicking phantom
formulation to produce phantoms with an absorption coefficient of
$0.05\text{\,}{\mathrm{cm}}^{-1}$ or $0.1\text{\,}{\mathrm{cm}}^{-1}$ at
$564\text{\,}\mathrm{nm}$ (the peak of the nigrosin spectrum).

####  2.1.2 Photoacoustic tomography (PAT)

For photoacoustic imaging, a small animal imaging system (MultiSpectral
Optoacoustic Tomography (MSOT) inVision 256-TF, iThera Medical) was used.
Briefly, a tunable optical parametric oscillator (OPO) pumped by an Nd:YAG
laser provides excitation pulses with a duration of $9\text{\,}\mathrm{ns}$ at
wavelengths from $660\text{\,}\mathrm{nm}$ to $1200\text{\,}\mathrm{nm}$ at a
repetition rate of $10\text{\,}\mathrm{Hz}$ with a wavelength tuning speed of
$10\text{\,}\mathrm{ms}$ and a peak pulse energy of $90\text{\,}\mathrm{mJ}$
at $720\text{\,}\mathrm{nm}$. Ten arms of a fibre bundle provide near-uniform
illumination over a disk extending approximately $8\text{\,}\mathrm{mm}$ along
the imaging chamber. Photoacoustic signals are detected using 256 toroidally
focused ultrasound transducers with a centre frequency of
$5\text{\,}\mathrm{MHz}$ (60% bandwidth), organized in a concave array of
$270\text{\,}\mathrm{\SIUnitSymbolDegree}$ angular coverage and a radius of
curvature of $4\text{\,}\mathrm{cm}$.

####  2.1.3 Detection of pO2

Two oxygen fluorescence quenching needle probes (NX-BF/O/E, Oxford Optronix)
were inserted into the flow circuit before and after the tissue-mimicking
phantom. A touch-screen monitor (OxyLite Pro, Oxford Optronix) displayed the
temperature and partial pressure of oxygen (pO2) real-time, and these data
were downloaded via an Arduino UNO and read in MATLAB. The relationship
between the oxygen saturation (sO2) and the partial pressure (pO2) in blood is
described by the characteristic sigmoid-shaped oxygen-hemoglobin dissociation
curve. A widely accepted fit to this curve is given by the Severinghaus
equation [31, 32]:

$\text{sO}_{2}(\%)=\frac{100}{\left[23400\times(\text{pO}_{2})^{3}+150\times(\text{pO}_{2})\right]^{-1}+1},$
(2)

which was used to convert our pO2 measurements into sO2.

####  2.1.4 Online flow spectrometry

Absorption spectra were recorded continuously via AvaSoft software using a
light source (Avalight-HAL-S-Mini, Avantes) and spectrometer (AvaSpec-
ULS2048-USB2-VA-50, Avantes) as the fluid passed through a flow cell
(170700-0.5-40, Hellma Analytics).

####  2.1.5 Offline spectrometry

To independently validate spectrophotometric measurements made in the flow
circuit, optical absorbance spectra were also recorded using a microplate
spectrometer (CLARIOstar, BMG LABTECH). Fluid samples were measured in a
96-well plate (Corning Costar).

###  2.2 Flow system characterization

####  2.2.1 Tubing assessment

The optimum tubing was determined by comparing photoacoustic images obtained
of agar phantoms containing various tube types filled with a
$25\text{\,}\mathrm{\SIUnitSymbolMicro mol}$ ICG solution (Sigma-Aldrich
I2633). The tubes were labelled according to their nominal inner and outer
diameters (I.D / O.D.) in $\text{\,}\mathrm{\SIUnitSymbolMicro m}$, but were
also made of different materials: polypropene 2660/2800 (Alliance Online
PSTS0007); THV500 2800/3150, 500/600 (Paradigm Optics); silicone 1570/2410,
630/1190, 300/630 (VWR 228-0256, 228-0254, 228-0253); PVC 1500/2100 (VWR
228-3857); PMMA 667/1000, 432/865, 375/500 (Paradigm Optics); polythene
580/960 (Smiths Medical 800/100/200 12665497).

####  2.2.2 Dye dilution series

A dilution series was used to demonstrate the utility of the flow circuit for
injecting different fluids into the closed PAT system, allowing PAT to be
performed concurrently with online spectrometry. Methylene blue (MB, 50484,
Fluka) and indocyanine green (ICG, I2633, Sigma-Aldrich) solutions were
prepared by diluting concentrations of $500\text{\,}\mathrm{\SIUnitSymbolMicro
M}$ and $100\text{\,}\mathrm{\SIUnitSymbolMicro M}$ (for MB and ICG
respectively) in deionized water. First, four spectra were measured and
averaged for samples of each concentration placed in a 96-well plate in the
CLARIOstar spectrophotometer. The concentrations were then flushed
individually through the flow system, starting with deionized water and then
sequentially with increasing concentration. Ten online spectra were recorded
over a range of $333\text{\,}\mathrm{nm}$ to $1100\text{\,}\mathrm{nm}$, and
ten single-slice PA images (no pulse-to-pulse averaging) were acquired for 17
wavelengths ($660\text{\,}\mathrm{nm}$, $664\text{\,}\mathrm{nm}$,
$680\text{\,}\mathrm{nm}$, $684\text{\,}\mathrm{nm}$,
$694\text{\,}\mathrm{nm}$, $700\text{\,}\mathrm{nm}$,
$708\text{\,}\mathrm{nm}$, $715\text{\,}\mathrm{nm}$,
$730\text{\,}\mathrm{nm}$, $735\text{\,}\mathrm{nm}$,
$760\text{\,}\mathrm{nm}$, $770\text{\,}\mathrm{nm}$,
$775\text{\,}\mathrm{nm}$, $779\text{\,}\mathrm{nm}$,
$800\text{\,}\mathrm{nm}$, $850\text{\,}\mathrm{nm}$,
$950\text{\,}\mathrm{nm}$), taking the mean of the ten single-slice images for
each concentration. The circuit was flushed with water between each
concentration.

####  2.2.3 Dynamic concentration change

A dynamic concentration change was used to illustrate the possibility for
real-time spectroscopic and photoacoustic measurements. Continuous acquisition
of online spectra and single-slice photoacoustic images was commenced once
deionized water was circulating within the flow system. After a certain time,
a high concentration dye solution (either
$500\text{\,}\mathrm{\SIUnitSymbolMicro M}$ MB or
$100\text{\,}\mathrm{\SIUnitSymbolMicro M}$ ICG) was injected using the
syringe pump at
($100\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{L}\mathrm{/}\min$).

The amount (moles) of the dye $Q(t)$ in the circuit at any given time $t$ can
be modelled using a first order differential equation. In Eq. (3), the rate of
change $dQ(t)/dt$ is equal to the difference between the inflow and outflow
amounts, expressed in terms of the flow rate $F$ induced by the syringe pump
($100\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{L}\mathrm{/}\min$), the
concentration $c$ of dye injected ($500\text{\,}\mathrm{\SIUnitSymbolMicro M}$
or $100\text{\,}\mathrm{\SIUnitSymbolMicro M}$) and the volume $V$
($5\text{\,}\mathrm{mL}$) of circulating fluid:

$\frac{dQ(t)}{dt}=F\left[c-\frac{Q(t)}{V}\right].$ (3)

The initial condition is given by Eq. (4), giving rise to Eq. (5), which is
the solution for the dye concentration $Q/V$ as a function of time. In the
limit $t\rightarrow\infty$, $Q(t)/V\rightarrow c$, as expected.

$Q(0)=Q_{0}=0;$ (4)

$\frac{Q(t)}{V}=c\left[1-e^{-\frac{F}{V}t}\right].$ (5)

###  2.3 Blood oxygenation measurements

Mouse blood, predominantly made up of strains B6 and 129SvEv, was collected
post mortem from the animal facility at the Cancer Research UK Cambridge
Institute, complying with the UK Animals (Scientific Procedures) Act 1986.
Each mouse contributed about $1\text{\,}\mathrm{mL}$ to a pool of about
$5\text{\,}\mathrm{mL}$ blood required for the flow circuit. Chemicals were
added to preserve, oxygenate and deoxygenate the blood. Prior to the
experiment, ethylenediaminetetraacetic acid (EDTA) anti-coagulant (9002-07-7,
Sigma-Aldrich) was added to the fresh blood, which was kept in a refrigerator
under $4\text{\,}\mathrm{\SIUnitSymbolCelsius}$ for no more than 72 hours.
During the experiment, the blood was returned to room temperature, and the
oxygenation was controlled chemically by adding $0.2\text{\,}\mathrm{\char
37\relax}$ v/v hydrogen peroxide, H2O2 (7722-84-1, Sigma-Aldrich), for
oxygenation [33] and sodium hydrosulfite (7775-14-6, ACROS Organics) for
deoxygenation [34]. The sodium hydrosulfite was dissolved in phosphate-
buffered saline (PBS) $\sim 0.03\text{\,}\mathrm{\char 37\relax}$ w/v for
injection into the flow circuit.

###  2.4 Spectral and statistical analysis

PAT contrast is provided by optical absorbers within the field of illumination
and detection. In general, PAT does not have sufficient resolution to
visualize individual molecules and therefore each image pixel (or voxel)
corresponds to more than one optical absorber, and a “spectral unmixing”
approach is required in order to extract the individual spectral components.
The linear mixture model [35] assumes that the measured spectrum is a linear
combination of distinct spectra:

$\textbf{x}=\textbf{SC}+n,$ (6)

where x is the $M$-wavelength $\times$ $N$-pixels measurement matrix, S is the
$M\times K$ matrix of $K$ spectra (“end-members”), C is the unknown $K\times
N$ matrix of end-member abundances (concentrations), and $n$ is measurement
noise. Neglecting $n$, an estimation Ĉ of the absorber concentrations can be
calculated from Eq. (6) by solving the following least-squares problem [11,
36]:

$\textbf{\^{C}}=\underset{C}{\arg\min}\left\lVert\textbf{SC}-\textbf{x}\right\rVert_{2}^{2},$
(7)

with the solution

$\textbf{\^{C}}=\textbf{xS}^{+},$ (8)

where $\textbf{S}^{+}$ is the pseudoinverse of S. The performance of spectral
unmixing carried out on acquired PA images was compared using three different
spectra: those from the literature [7, 8, 9], and those experimentally
measured in this study from either the online spectrometer cell or the offline
plate reader. In addition, the online flow spectra were used to provide an
independent ground truth for sO2; this entailed unmixing of the online flow
spectra with the spectra for oxy- and deoxy- hemoglobin also measured live
during the experiments.

To correct for the spectrally varying fluence (“spectral coloring”) when using
optically absorbing phantoms, the image intensities in the tube region were
divided by the nigrosin absorption spectrum. Specifically, the nigrosin
absorption spectrum was normalised to the known absorption coefficient (at the
$564\text{\,}\mathrm{nm}$ peak) and re-sampled at the 17 experimental
wavelengths; the multispectral PA data were then divided by a wavelength-
specific factor calculated from exponential decay over the
$8.5\text{\,}\mathrm{mm}$ distance (the background agar material impregnated
with absorbing nigrosin) from the outer edge of the phantom to the tube wall.
All data and source code used in this publication is available on
https://doi.org/10.17863/CAM.40365 and GitHub (https://github.com/9xg/flow-
phantom).

##  3 Results

###  3.1 Flow system characterization

Eleven different tubes were assessed for their photoacoustic imaging
suitability (Fig. 2). Tubes made from PMMA (I.D./O.D. 375/500, 432/865,
667/1000) showed the lowest signal-to-background (SBR) ratios ranging from
0.15 to 0.37. The other tubes materials show a clear trend of increasing SBR
with larger inner diameters. Silicone tubes present high SBRs, even for low
I.D./O.D. ratios (300/600, 630/1190, 1570/2410). Three PVC tubes showed very
similar performance (I.D./O.D. 1500/2100, 1570/2410, and 2660/2800). The PVC
tube with 1500/2100 presented a uniform, circular appearance in images with
high SBR. It was therefore selected for the remaining experiments.

Figure 2: Assessment of the effect of tube type on the quality of
photoacoustic images. (a) PAT images of four tubes with inner and outer
diameters (I.D.-O.D.) in $\text{\,}\mathrm{\SIUnitSymbolMicro m}$ of 630-1190
(silicone), 667-1000 (PMMA), 1500-2100 (PVC) and 1570-2410 (silicone), filled
with a solution of ICG and embedded within a scattering agar cylinder (not
shown). (b) Mean photoacoustic signal intensity inside the tube relative to
the signal outside (signal-to-background ratio, SBR) for eleven different
tubes. The tube with I.D.-O.D.
$1500\text{\,}$-$2100\text{\,}\mathrm{\SIUnitSymbolMicro m}$ was selected due
to its high SBR, flexibility and low cost. The four tubes illustrated in (a)
are marked with symbols.

Dilutions of MB and ICG were tested inside the flow phantom under closed
conditions i.e. with a steady concentration flowing in the circuit. By
calculating the mean PA signal intensities over a manually segmented tube
cross-sectional region of interest (ROI) for 17 different wavelengths it was
possible to compare the PA spectra with those acquired during online flow
spectrometry (Fig. 3 (a-d)). Unexpectedly, the two independent spectral
measurements show poor agreement in terms of relative intensities and spectral
shape for the different dye concentrations. In particular, it is notable that
the PA signal intensities (Fig. 3b) show a more pronounced shift in spectral
peak with increasing ICG concentration compared with the online spectrometer
data (Fig. 3d); this may be a result of spectral coloring where light
attenuation at the $800\text{\,}\mathrm{nm}$ peak suppresses the PA signal
intensity deeper within the tube, leading to an overall reduction in mean
intensity. For MB the discrepancies are difficult to discern seeing as the
prominent spectral features occur below the PA imaging range, but it is
interesting that the relative PA intensities are about half those for ICG
suggesting that MB has poorer PA signal generation efficiency. Overall, these
spectral inconsistencies raise a question about which are the most suitable
endmember spectra for unmixing PAT images of MB and ICG.

Figure 3: Comparison of dye spectra measured using PAT and online flow
spectrometry under closed conditions. (a,b) PA signal intensities (with
smoothed spline) and (c,d) online flow spectrometer absorbance values measured
for different concentrations of MB (a,c), and ICG (b,d). Grey box indicates
wavelengths outside the PAT spectral range.

The spectra measured with the online spectrometer and offline in the plate
reader were used to perform concentration-specific unmixing of the mean ROI
intensities, and these unmixed intensities were compared with those calculated
using the literature spectra incorporated within the PA analysis software.
These literature spectra are shown in Fig. 4 (a-b) along with online and
offline spectra for example concentrations (Fig. 4 (c-f)).

Figure 4: Comparison of endmember spectra used for spectral unmixing.
Literature spectra (a,b), spectra acquired using the online flow spectrometer
(c,d) and the offline flow spectrometer (e,f) for MB (a,c,e) and ICG (b,d,f).

Figure 5 (a-b) shows a linear relationship between dye concentration and
spectrally unmixed PA intensities up to
$100\text{\,}\mathrm{\SIUnitSymbolMicro mol}$ for MB, and
$40\text{\,}\mathrm{\SIUnitSymbolMicro mol}$ for ICG, respectively. Divergence
between the two types of unmixing is likely due to concentration dependent
changes in the spectra of these dyes, which are particularly prominent for
ICG.

Figure 5: Spectrally unmixed PA signal intensities for a range of dye
concentrations (a) MB and (b) ICG concentrations obtained by unmixing with
literature spectra and those measured using the online and offline
spectrometers. The linear model (LM) was used to calculate concentrations of
MB (c) and ICG (d) during evolution of dye concentrations within the flow
circuit due to continual dye injection. The starting points were forced to
zero for ease of comparison.

Having explored the application of the flow circuit under static concentration
values, we then examined the response of the PAT instrument to dynamic changes
in dye concentrations for MB and ICG (Fig. 5 (c-d)). PA images were unmixed
using the literature spectra and those recorded live using the online
spectrometer, and then converted to absolute concentration values using the
linear models calculated and plotted in (a-b). Neglecting the initial lag
phase (which could not be experimentally determined and therefore was not
incorporated into the model), the rates of concentration change match
reasonably well with those predicted by the model (Eq. 5). However, there is a
notable discrepancy between the dynamics calculated using the literature and
online spectra for ICG unmixing; if the lag phase were corrected for, the
online unmixing would be the closer match to the model and this is again
likely to be a consequence of the concentration-dependent change in spectral
shape, which is incorporated in the online spectra but not the single
literature spectrum.

###  3.2 Blood oxygenation measurements

To demonstrate the closed nature of the blood flow circuit, blood was fully
oxygenated by adding approximately $15\text{\,}\mathrm{\SIUnitSymbolMicro L}$
$0.2\text{\,}\mathrm{\char 37\relax}$ v/v H2O2 to $8\text{\,}\mathrm{mL}$
mouse blood and then injecting this into the flow system and circulating it
for several minutes. Figure 6 shows the blood oxygen saturation (sO2) measured
within the flow circuit by: the pO2 probe; the online spectrometer and PAT.
The sO2 was calculated from the pO2 data using Equation 2 and from the online
spectrometer and PAT data through spectral unmixing. To enable comparison
between the methods independently, unmixing was performed using the spectra
inherent to each system: for the online spectrometer, experimentally measured
spectra were used as endmembers, while for PAT data, the literature spectra
were used. As hoped, all three sO2 measurements show consistent values over
time. However, the sO2 calculated from the unmixed PA images under-reads by
about $13\text{\,}\mathrm{\char 37\relax}$ compared to the ground truth sO2
calculated from the pO2 probe and the online flow spectrometer.

Figure 6: Assessment of blood oxygenation within the flow circuit under closed
conditions. Blood oxygen saturation (sO2) was calculated using the three
independent methods: the pO2 probe measurements (using Equation 2); unmixing
of the spectra measured using the online flow spectrometer (with
experimentally measured spectra as end-members); unmixing of the mean pixel
intensities in the PA images (with literature spectra as end-members) acquired
while blood circulated in the flow system for 200 seconds.

In addition to investigating static blood oxygenation within the circuit, the
full range of blood sO2 values ($100\text{\,}\mathrm{\char 37\relax}$ to
$0\text{\,}\mathrm{\char 37\relax}$) were explored by injecting sodium
hydrosulfite into the flow circuit in order to gradually deoxygenate the
blood. The blood sO2 values cover the expected dynamic range when calculated
from both the pO2 probe and the online spectrometer (Fig. 7a); however, PA
signal intensities show markedly different behaviors when unmixed with spectra
(Fig. 7b) from the literature or from the online flow spectrometer. The sO2
values calculated by unmixing with the experimentally measured spectra (red
line, Fig. 7a) show good agreement with the spectrometer-derived ground truth
values (yellow dotted line, Fig. 7a), except for somewhat over-reading the
oxygenation for sO2 values below about $20\text{\,}\mathrm{\char 37\relax}$.
They also show a dynamic range comparable to that obtained from the pO2 probe.
However, the sO2 values calculated by unmixing with the literature spectra
surprisingly exhibit a dramatically reduced dynamic range
($86\text{\,}\mathrm{\char 37\relax}$ to $58\text{\,}\mathrm{\char
37\relax}$); this may be due to the unmixing spectra and/or spectral coloring.

Figure 7: Dynamic deoxygenation of the circulating blood. (a) Change in blood
oxygen saturation (sO2) measured using the pO2 probe, online spectrometer and
PAT system while injecting 3% w/v sodium hydrosulfite in PBS over a period of
nine minutes. (b) Literature spectra for oxygenated (HbO2) and deoxygenated
(Hb) hemoglobin used for spectral unmixing [7] compared to spectra obtained
from the online spectrometer at the start and end of the dynamic
deoxygenation.

Two further investigations explored the potential impact of spectral
colouring. First, sO2 measurements across the tube diameter were calculated
(Fig. 8a). Unlike the results in Fig. 6 and Fig. 7 which were calculated from
the mean of the ROI defined on the tube cross-sectional area, the results in
Fig. 8a entailed radial unmixing using the mean of pixels found in a circle
around the perimeter of the ROI, and repeating this for circle diameters
reducing in one pixel ($75\text{\,}\mathrm{\SIUnitSymbolMicro m}$) increments.
For the fully oxygenated blood the sO2 values are consistent across the tube
diameter, but as the blood becomes deoxygenated there is a clear trend towards
sO2 over-reading at the centre of the tube. This is a clear demonstration of
spectral coloring. A second illustration of spectral coloring is shown in Fig.
8b where experiments were repeated with absorbing nigrosin dye incoporated
into the background of the tissue mimicking phantom surrounding the flow
circuit tube. Increasing the background absorption causes increasing
inaccuracy of the sO2 relative to the ground truth since wavelengths of light
are preferentially absorbed by the nigrosin dye, leading to over-estimation of
the contribution of oxyhemoglobin to the PA signal. A simple light fluence
correction, implemented by dividing the images by the known nigrosin spectrum
(section 2.4), restored consistent sO2 measurements irrespective of the
background absorption, providing evidence that this spatial phenomenon is
indeed related to spectral coloring.

Figure 8: Investigation of spectral coloring. (a) Spatial profiles of blood
oxygen saturation (sO2) as a function of radial position within the flow
circuit tube by PA unmixing with spectra measured using the online
spectrometer during dynamic blood oxygenation (Fig. 7). (b) Impact of
correction for the absorption spectrum of the nigrosin dye included in the
background of the tissue mimicking phantom as a function of dye concentration.
For absorption coefficients of $>$0\text{\,}{\mathrm{cm}}^{-1}$$ the sO2 (grey
bars) deviated from the ground truth value (dotted line); accurate values were
restored by recalculating sO2 (orange bars) after dividing the photoacoustic
images by the relevant nigrosin spectrum as described in the text.

##  4 Discussion

The application of PAT in assessment of blood hemoglobin concentration and
oxygen saturation is now widespread, but estimation of the oxy- and deoxy-
hemoglobin chromophore concentrations is affected by a variety of tissue- and
instrumentation-dependent factors. We have developed and applied a low-cost
flow phantom system with online monitoring of optical absorption spectra and
partial pressure of oxygen in order to facilitate testing of PAT systems.

We first selected an appropriate tubing to maximise PAT image quality and
tested the circuit using two optically absorbing dyes. For both dyes, MB and
ICG, the absorption spectrum exhibits a double peak where the relative peak
intensities shift with increasing concentration. Unmixing with a single
“average” spectrum is therefore unlikely to yield the same result as unmixing
with a concentration-specific (online/offline) spectrum; indeed, we observed
different and slightly non-linear relationships between dye concentration and
spectrally unmixed intensities (which are assumed to be proportional to
concentration) for the two types of unmixing. For ICG, the values from
unmixing with the literature spectrum under-read those from unmixing with the
online/offline spectra (Fig. 5b,d); for MB the unmixing results were
indistinguishable (Fig. 5a,c) since the concentration dependent effects occur
outside of the wavelength range of the PA system under test. It is possible
that erroneous concentration measurements may also be caused by spectral
coloring, even though only a single component is used in the spectral
unmixing; two component unmixing (Hb and HbO2) is more relevant for
calculation of sO2.

We proceeded to evaluate the estimation of sO2 under both closed circuit and
dynamic conditions. The results highlight two major factors that contribute to
inaccuracies in spectral unmixing, and therefore in sO2 measurements. The
first is the choice of spectra for unmixing the photoacoustic images; the
second is the effect of spectral coloring. The online spectrometer allowed
spectra to be recorded during PAT image acquisition and used as alternatives
to published spectra for unmixing the images. Moreover, the ground truth sO2
values provided by the spectrometer and the pO2 probes are designed to
facilitate exploration of phenomena such as spectral coloring.

Images of blood flowing through the circuit entailed unmixing with two
components: oxy- and deoxy- hemoglobin. The resulting sO2 calculations
(HbO2/THb) were consistent over time but were dependent on the hemoglobin
spectra used for unmixing, the blood oxygenation saturation, and also the
position within the tube. The spectra for oxy- and deoxy- hemoglobin measured
online were considerably different from those available in the literature that
are widely used for spectral unmixing both in vitro and in vivo. Using our
online flow spectrometer spectra improved the accuracy of the sO2 calculations
compared to those made using the default spectra, suggesting that the online
spectra were more representative of the actual optical absorption presented to
the PAT system.

Over-reading of the lowest sO2 values was also observed and attributed to
spectral coloring: the high absorption by deoxygenated blood at wavelengths
below about $750\text{\,}\mathrm{nm}$ results in low light intensities at
these wavelengths in the centre of the tube and therefore the PA intensities
(proportional to both absorption and light fluence) varied with wavelength in
a way more closely matching with oxyhemoglobin, leading to overestimation of
the oxyhemoglobin concentration. Admittedly, these sO2 values are well below
the physiological range ($70\text{\,}\mathrm{\char
37\relax}100\text{\,}\mathrm{\char 37\relax}$), nonetheless it is important to
explore the full dynamic range in order to establish the accuracy limits of PA
sO2 measurements. A simple example of correcting the effects of spectal
coloring introduced by the tissue mimicking phantom surrounding the tube was
achieved through dividing by the background nigrosin spectrum, which restored
accurate sO2 measurements.

One limitation of the study is the discrepancy between our measured hemoglobin
spectra and those widely used in the literature. This could be attributed to
optical scattering due to the blood cells within the measurement cuvette, as
literature studies normally measure the haemoglobin molecule directly; another
factor could be the strong attenuation of the blood, necessitating use of a
neutral density filter in order to prevent signal saturation when recording
the reference spectrum (taken with PBS in the tube), but then removal of the
filter for the blood measurements. Although the filter was taken into account
in the calculation of abosrption spectra, the effect of blood scattering is
difficult to compensate for; future experiments should therefore calibrate the
measured spectra by extracting blood samples from the circuit at different sO2
values and validating the sO2 using a blood gas analyser. It would also be
interesting to investigate oxygenation states of haemoglobin extracted from
blood cells thereby obviating optical scattering effects. A further limitation
arises from the pO2 measurements, as the dynamics of the deoxygenation study
did not directly mirror the sO2 values obtained in the spectrometer. This may
be due to an unrepresentative conversion from pO2 to sO2: the Severinghaus
equation is derived from human, not mouse, blood data and also assumes certain
values for parameters such as pH and temperature. Alternative conversions such
as the Kelman equation were also explored, and measured values for pH and
temperature were incorporated, but this did not significantly alter the trend
in sO2 over time. It is possible that bubbles or insufficient contact with the
blood in the circuit corrupted the pO2 readings, and therefore future versions
of the circuit will integrate truly in-flow pO2 probes and a membrane
oxygenator to avoid bubbles.

The application of the presented phantom opens a range of opportunities for
future studies of tissue- and instrument-dependent correction factors in PAT.
In particular, future work will attempt to model and account for the spectral
distortion across the tube diameter observed during blood oxygenation studies.
The development of corrections incorporating multiple choromphores inside and
outside the tube both ex vivo and in vivo is more complex and continues to be
investigated, for example using Monte Carlo simulations and model-based
iterative minimisation [16, 17, 21, 22, 23, 19]. Moreover, the number and
choice of wavelengths used for spectral unmixing is an important consideration
[13] to be investigated in future work.

##  5 Conclusion

In summary, we have developed a low-cost flow phantom that includes an online
spectrometer and partial pressure of oxygen probe to facilitate detailed
validation of PAT measurements of blood oxygen saturation, sO2. We found that
it is important to correctly identify the absorption spectra to be used for
unmixing photoacoustic images in order to accurately determine even relative
absorber concentrations. We also found that the co-dependence of photoacoustic
signal intensity on light fluence and absorption leads to a major challenge in
accounting for spectral coloring, which can lead to substantial
underestimation of sO2. Further calibration and automation of the circuit will
enable additional ex vivo studies requiring careful control and knowledge of
blood sO2, as well as opening the possibility of exploring additional blood
parameters such as pH, and also absorptive and other properties of different
fluids. Detailed understanding of the photoacoustic signal origins ex vivo
remains essential for proper interpretation of photoacoustic measurements made
in vivo.

### Disclosures

SEB has received research support from iThera Medical GmbH and PreXion Inc.,
vendors of photoacoustic imaging instruments.

###### Acknowledgements.

The authors would like to thank Ayaka Shinozaki for her assistance with data
collection, James Joseph for his contributions to the experimental design, and
Michael Schneider for his helpful input regarding data interpretation. This
work was supported by Cancer Research UK (C47594/A16267, C14303/A17197) and
the EPSRC-CRUK Cancer Imaging Centre in Cambridge and Manchester (C197/A16465
and C8742/A18097). We would like to thank the CRUK CI Core Facilities for
their support of this work, in particular the Imaging Core and the Biological
Resource Unit.

## References

  * [1] M. Martinho Costa, A. Shah, I. Rivens, et al., “Quantitative photoacoustic imaging study of tumours in vivo: Baseline variations in quantitative measurements,” Photoacoustics 13, 53–65 (2019).
  * [2] I. Quiros-Gonzalez, M. R. Tomaszewski, S. J. Aitken, et al., “Optoacoustics delineates murine breast cancer models displaying angiogenesis and vascular mimicry,” British Journal of Cancer 118, 1098–1106 (2018).
  * [3] J. Yang, G. Zhang, Q. Li, et al., “Photoacoustic imaging for the evaluation of early tumor response to antivascular treatment,” Quantitative Imaging in Medicine and Surgery 9(2), 160–170 (2019).
  * [4] M. J. Waldner, F. Knieling, C. Egger, et al., “Multispectral Optoacoustic Tomography in Crohn’s Disease: Noninvasive Imaging of Disease Activity,” Gastroenterology 151(2), 238–240 (2016).
  * [5] N. Bhutiani, W. E. Grizzle, S. Galandiuk, et al., “Noninvasive Imaging of Colitis Using Multispectral Optoacoustic Tomography,” Journal of Nuclear Medicine 58(6), 1009–1012 (2017).
  * [6] H. F. Zhang, K. Maslov, M. Sivaramakrishnan, et al., “Imaging of hemoglobin oxygen saturation variations in single vessels in vivo using photoacoustic microscopy,” Applied Physics Letters 90(5), 053901 (2007).
  * [7] S. Prahl, Optical Absorption of Hemoglobin (1999).  https://omlc.org/spectra/hemoglobin/.
  * [8] S. Prahl, Methylene Blue Spectra (2017).  https://omlc.org/spectra/mb/.
  * [9] M. Landsman, G. Kwant, G. Mook, et al., “Light-absorbing properties, stability, and spectral stabilization of indocyanine green,” Journal of Applied Physiology 40(4), 575–583 (1976).
  * [10] K. Maslov, H. F. Zhang, and L. V. Wang, “Effects of wavelength-dependent fluence attenuation on the noninvasive photoacoustic imaging of hemoglobin oxygen saturation in subcutaneous vasculature in vivo,” Inverse Problems 23, S113–S122 (2007).
  * [11] B. Cox, J. G. Laufer, S. R. Arridge, et al., “Quantitative spectroscopic photoacoustic imaging: a review.,” Journal of Biomedical Optics 17(6), 061202 (2012).
  * [12] R. O. Esenaliev, I. V. Larina, K. V. Larin, et al., “Optoacoustic technique for noninvasive monitoring of blood oxygenation: a feasibility study,” Applied Optics 41(22), 4722 (2002).
  * [13] R. Hochuli, P. C. Beard, and B. Cox, “Effect of wavelength selection on the accuracy of blood oxygen saturation estimates obtained from photoacoustic images,” in Proc. of SPIE, Photons Plus Ultrasound: Imaging and Sensing, A. A. Oraevsky and L. V. Wang, Eds., 9323, 93231V–8 (2015).
  * [14] T. Mitcham, H. Taghavi, J. Long, et al., “Photoacoustic-based sO2 estimation through excised bovine prostate tissue with interstitial light delivery.,” Photoacoustics 7, 47–56 (2017).
  * [15] J. Laufer, C. Elwell, D. Delpy, et al., “In vitro measurements of absolute blood oxygen saturation using pulsed near-infrared photoacoustic spectroscopy: accuracy and resolution.,” Physics in Medicine and Biology 50, 4409–4428 (2005).
  * [16] B. Cox, S. Arridge, K. Köstli, et al., “Two-dimensional quantitative photoacoustic image reconstruction of absorption distributions in scattering media by use of a simple iterative method,” Applied Optics 45(8), 4 – 6 (2006).
  * [17] Z. Yuan and H. Jiang, “Quantitative photoacoustic tomography: Recovery of optical absorption coefficient maps of heterogeneous media,” Applied Physics Letters 88(23), 231101 (2006).
  * [18] J. Laufer, D. Delpy, C. Elwell, et al., “Quantitative spatially resolved measurement of tissue chromophore concentrations using photoacoustic spectroscopy: application to the measurement of blood oxygenation and haemoglobin concentration,” Physics in Medicine and Biology 52, 141–168 (2007).
  * [19] J. Laufer, B. Cox, E. Zhang, et al., “Quantitative determination of chromophore concentrations from 2D photoacoustic images using a nonlinear model-based inversion scheme.,” Applied Optics 49, 1219–33 (2010).
  * [20] Y. Liu, H. Jiang, and Z. Yuan, “Two schemes for quantitative photoacoustic tomography based on Monte Carlo simulation,” Medical Physics 43(7), 3987–3997 (2016).
  * [21] L. Yao, Y. Sun, and H. Jiang, “Quantitative photoacoustic tomography based on the radiative transfer equation.,” Optics Letters 34, 1765–7 (2009).
  * [22] F. M. Brochu, J. Brunker, J. Joseph, et al., “Towards Quantitative Evaluation of Tissue Absorption Coefficients Using Light Fluence Correction in Optoacoustic Tomography,” IEEE Transactions on Medical Imaging 36, 322–331 (2017).
  * [23] B. T. Cox, S. R. Arridge, and P. C. Beard, “Estimating chromophore distributions from multiwavelength photoacoustic images,” Journal of Optical Society of America A 26(2), 443–455 (2009).
  * [24] S. Tzoumas, A. Nunes, I. Olefir, et al., “Eigenspectra optoacoustic tomography achieves quantitative blood oxygenation imaging deep in tissues,” Nature Communications 7(May), 12121 (2016).
  * [25] W. C. Vogt, C. Jia, K. A. Wear, et al., “Biologically relevant photoacoustic imaging phantoms with tunable optical and acoustic properties,” Journal of Biomedical Optics 21(10), 101405 (2016).
  * [26] E. Maneas, W. Xia, O. Ogunlade, et al., “Gel wax-based tissue-mimicking phantoms for multispectral photoacoustic imaging,” Biomedical Optics Express 9(3), 1151 (2018).
  * [27] W. C. Vogt, X. Zhou, R. Andriani, et al., “Photoacoustic oximetry imaging performance evaluation using dynamic blood flow phantoms with tunable oxygen saturation,” Biomedical Optics Express 10(2), 449 (2019).
  * [28] S. N. Hennen, W. Xing, Y. B. Shui, et al., “Photoacoustic tomography imaging and estimation of oxygen saturation of hemoglobin in ocular tissue of rabbits,” Experimental Eye Research 138, 153–158 (2015).
  * [29] Z. Chen, S. Yang, and D. Xing, “In vivo detection of hemoglobin oxygen saturation and carboxyhemoglobin saturation with multiwavelength photoacoustic microscopy,” Optics Letters 37(16), 3414 (2012).
  * [30] Y. Wang, S. Hu, K. Maslov, et al., “In vivo integrated photoacoustic and confocal microscopy of hemoglobin oxygen saturation and oxygen partial pressure,” Optics Letters 36(7), 1029 (2011).
  * [31] J. W. Severinghaus, “Simple, accurate equations for human blood O2 dissociation computations .,” Journal of Applied Physiology: Re 46(3), 599–602 (1979).
  * [32] J.-A. Collins, A. Rudenski, J. Gibson, et al., “Relating oxygen partial pressure, saturation and content: the haemoglobin-oxygen dissociation curve,” Breathe 11, 194–201 (2015).
  * [33] D. White and P. Teasdale, “The oxygenation of blood by hydrogen peroxide,” British Journal of Anaesthesia 36, 528–529 (1964).
  * [34] K. Briely-Sabo and A. Bjornerud, “Accurate de-oxygenation of ex vivo whole blood using sodium dithionite,” Proc. Intl. Sot. Mag. Reson. Med 117(1985), 2025 (2000).
  * [35] S. Tzoumas, N. Deliolanis, S. Morscher, et al., “Unmixing molecular agents from absorbing tissue in multispectral optoacoustic tomography,” IEEE Transactions on Medical Imaging 33(1), 48–60 (2014).
  * [36] L. Ding, X. L. Dean-Ben, N. C. Burton, et al., “Constrained Inversion and Spectral Unmixing in Multispectral Optoacoustic Tomography,” IEEE Transactions on Medical Imaging 36, 1676–1685 (2017).

Marcel Gehrung completed his MSc degree at Eberhard Karls University Tuebingen
in 2017 and is presently completing his PhD in Medical Science at the
University of Cambridge focused on the application of advanced image analysis
methods for cancer diagnosis.

Sarah Bohndiek completed her BA degree at the University of Cambridge in 2005
and her PhD in Radiation Physics at University College London in 2008,
specialising in X-ray diffraction. After postdoctoral fellowships in both the
UK and USA, she started the VISIONLab in Cambridge in 2013, which develops and
applies new imaging biomarkers to shed light on the tumour microenvironment.

Joanna Brunker received her MSc degree in natural sciences (2009) and her PhD
(2013) in Medical Physics and Biomedical Engineering from the University
College London, UK. Following her postdoctoral research fellowships at UCL and
at the Cancer Research UK Cambridge Institute, she is now (July 2019) starting
a research group at the UCL Wellcome / EPSRC Centre for Interventional and
Surgical Sciences (WEISS) where she will focus on translation of photoacoustic
and other imaging technologies into clinical practice.

###### List of Figures

  1. 1 Overview of the flow system. (1) Injection site for introducing oxygenated blood (or other fluids) into the flow system, and for subsequently deoxygenating the blood using sodium hydrosulfite delivered via the syringe driver (MKCB2159V, Harvard); (2) spectra are recorded using a light source (Avalight-HAL-S-Mini, Avantes) and spectrometer (AvaSpec-ULS2048-USB2-VA-50, Avantes) as the blood passes through a flow cell (170700-0.5-40, Hellma Analytics); (3) needle probes (NX-BF/O/E, Oxford Optronix) measure the temperature and partial pressure of oxygen (pO2) before (3a) and after (3b) the blood passes through an agar phantom immersed in the photoacoustic imaging system (MSOT inVision 256-TF, iThera Medical) (4); (5) a touch-screen monitor (OxyLite Pro, Oxford Optronix) displays temperature and oxygen data; (6) these data are downloaded via an Arduino UNO and read in MATLAB on a laptop, which also records the spectrometer readings via AvaSoft software; (7) a peristaltic pump (CTP100, Fisher Scientific) provides blood circulation.
  2. 2 Assessment of the effect of tube type on the quality of photoacoustic images. (a) PAT images of four tubes with inner and outer diameters (I.D.-O.D.) in $\text{\,}\mathrm{\SIUnitSymbolMicro m}$ of 630-1190 (silicone), 667-1000 (PMMA), 1500-2100 (PVC) and 1570-2410 (silicone), filled with a solution of ICG and embedded within a scattering agar cylinder (not shown). (b) Mean photoacoustic signal intensity inside the tube relative to the signal outside (signal-to-background ratio, SBR) for eleven different tubes. The tube with I.D.-O.D. $1500\text{\,}$-$2100\text{\,}\mathrm{\SIUnitSymbolMicro m}$ was selected due to its high SBR, flexibility and low cost. The four tubes illustrated in (a) are marked with symbols.
  3. 3 Comparison of dye spectra measured using PAT and online flow spectrometry under closed conditions. (a,b) PA signal intensities (with smoothed spline) and (c,d) online flow spectrometer absorbance values measured for different concentrations of MB (a,c), and ICG (b,d). Grey box indicates wavelengths outside the PAT spectral range.
  4. 4 Comparison of endmember spectra used for spectral unmixing. Literature spectra (a,b), spectra acquired using the online flow spectrometer (c,d) and the offline flow spectrometer (e,f) for MB (a,c,e) and ICG (b,d,f).
  5. 5 Spectrally unmixed PA signal intensities for a range of dye concentrations (a) MB and (b) ICG concentrations obtained by unmixing with literature spectra and those measured using the online and offline spectrometers. The linear model (LM) was used to calculate concentrations of MB (c) and ICG (d) during evolution of dye concentrations within the flow circuit due to continual dye injection. The starting points were forced to zero for ease of comparison.
  6. 6 Assessment of blood oxygenation within the flow circuit under closed conditions. Blood oxygen saturation (sO2) was calculated using the three independent methods: the pO2 probe measurements (using Equation 2); unmixing of the spectra measured using the online flow spectrometer (with experimentally measured spectra as end-members); unmixing of the mean pixel intensities in the PA images (with literature spectra as end-members) acquired while blood circulated in the flow system for 200 seconds.
  7. 7 Dynamic deoxygenation of the circulating blood. (a) Change in blood oxygen saturation (sO2) measured using the pO2 probe, online spectrometer and PAT system while injecting 3% w/v sodium hydrosulfite in PBS over a period of nine minutes. (b) Literature spectra for oxygenated (HbO2) and deoxygenated (Hb) hemoglobin used for spectral unmixing [7] compared to spectra obtained from the online spectrometer at the start and end of the dynamic deoxygenation.
  8. 8 Investigation of spectral coloring. (a) Spatial profiles of blood oxygen saturation (sO2) as a function of radial position within the flow circuit tube by PA unmixing with spectra measured using the online spectrometer during dynamic blood oxygenation (Fig. 7). (b) Impact of correction for the absorption spectrum of the nigrosin dye included in the background of the tissue mimicking phantom as a function of dye concentration. For absorption coefficients of $>$0\text{\,}{\mathrm{cm}}^{-1}$$ the sO2 (grey bars) deviated from the ground truth value (dotted line); accurate values were restored by recalculating sO2 (orange bars) after dividing the photoacoustic images by the relevant nigrosin spectrum as described in the text.


11institutetext: 1Institute of Artificial Intelligence, Beihang University,
Beijing, China
2The University of Hong Kong, Hong Kong, China
3Deepwise AI Lab, Beijing, China
4Hwa Mei Hospital, UCAS
5Ningbo Institute of Life and Health Industry, UCAS
6School of Artificial Intelligence, Xidian University, Xi’an, China
7School of Automation, Northwestern Polytechnical University, Xi’an, China

# Computer-aided Tuberculosis Diagnosis with Attribute Reasoning Assistance

Chengwei Pan 11 Gangming Zhao C. Pan and G. Zhao contributed equally to this
work.2233 Junjie Fang 4455 Baolian Qi 55 Jiaheng Liu 11
Chaowei Fang 66 Dingwen Zhang 77 Jinpeng Li${}^{\left(\textrm{{\char
0\relax}}\right)}$ Corresponding author. Email<EMAIL_ADDRESS>Yizhou
Yu 22 3 3

###### Abstract

Although deep learning algorithms have been intensively developed for
computer-aided tuberculosis diagnosis (CTD), they mainly depend on carefully
annotated datasets, leading to much time and resource consumption. Weakly
supervised learning (WSL), which leverages coarse-grained labels to accomplish
fine-grained tasks, has the potential to solve this problem. In this paper, we
first propose a new large-scale tuberculosis (TB) chest X-ray dataset, namely
tuberculosis chest X-ray attribute dataset (TBX-Att), and then establish an
attribute-assisted weakly supervised framework to classify and localize TB by
leveraging the attribute information to overcome the insufficiency of
supervision in WSL scenarios. Specifically, first, the TBX-Att dataset
contains 2000 X-ray images with seven kinds of attributes for TB relational
reasoning, which are annotated by experienced radiologists. It also includes
the public TBX11K dataset with 11200 X-ray images to facilitate weakly
supervised detection. Second, we exploit a multi-scale feature interaction
model for TB area classification and detection with attribute relational
reasoning. The proposed model is evaluated on the TBX-Att dataset and will
serve as a solid baseline for future research. The code and data will be
available at https://github.com/GangmingZhao/tb-attribute-weak-localization.

###### Keywords:

Computer-aided tuberculosis diagnosis Weakly supervised learning Attribute
reasoning

##  1 Introduction

Tuberculosis (TB) is one of the most serious thoracic diseases associated with
a high death rate. The early diagnosis and treatment is very important to
confront the threat of TB. Computer-aided tuberculosis diagnosis (CTD) has
been widely investigated to assist radiologists in diagnosing TB. Although
deep learning algorithms have achieved stupendous success in automatic disease
classification [24, 1] and localization [3, 11, 14] for chest X-rays, their
performance on CTD remains a barrier to clinical application. One obstacle
lies in the absence of high quality annotated datasets and the scarcity of
clinical features (attributes). In practice, such attributes have been
validated to be helpful to the performance and interpretability.

In addition, the main challenges of diagnosing TB with chest X-ray images
include the low visual contrast between lesions and other regions, and
distortions induced by other overlapping tissues. Sometimes it is difficult
for radiologists to recognize obscure diseases.For CTD, previous methods
mainly concentrated on disease classification [19, 2] , and several recent
works have taken a step forward to detect disease regions under weak/limited
supervisions. They can be grouped into two main categories: the first category
[24] resorts to convolutional neural networks (CNN) trained on the
classification task and output disease localization results through
calculating the category activation maps [29]; the second category [11, 14]
uses the multiple instance learning to directly yield categorical probability
maps that can be easily transformed to lesion positions. However, the
performance of these methods is still far from clinical usage. Effective WSL
algorithms remain the boundary to explore.

The mainstream pipeline of object detection is screening out potential
proposals followed by proposal classification [21]. By stacking piles of
convolutional layers, CNNs are very advantageous at extracting surrounding
contextual information, however, distant relationships are still hard to
capture by convolutions with small kernels. On the other hand, radiologists
use attribute information for diagnosis, whereas the pure data-driven CNNs can
not take this prior knowledge into account, leading to disconnections to the
clinical practice.

The main contributions of this paper can be summarized as follows.

  * •

We present a large-scale TBX-Att dataset for attribute-assisted X-ray
diagnosis for TB. To the best of our knowledge, this offers the first
playground to leverage the attribute information to help models detect TB
areas in a weakly-supervised manner.

  * •

We propose an effective method to fuse the attribute information and TB
information. The multi-scale attribute features are extracted for each TB
proposal under the guidance of a relational reasoning module and are further
refined with the attribute-guided knowledge. A novel feature fusion module is
devised to enhance the TB representation with its effective attribute prompt.

  * •

The proposed method improves object detection baselines [21, 12] by large
margins on TBX-Att, leading to a solid benchmark for weakly supervised TB
detection.

##  2 Related Work

###  2.1 Object Detection

Object detection is a widely-studied topic in both natural and medical images.
It aims at localizing object instances of interest such as faces, pedestrians
and disease lesions. The most famous kind of deep learning approaches for
object detection is the R-CNN [7] family. The primitive R-CNN extracts
proposals through selective search [23], and then predicts object bounding
boxes and categories from convolution features of these proposals. Fast R-CNN
[6] adopts a shared backbone network to extract proposal features via RoI
pooling. Faster R-CNN [21] automatically produces object proposals from top-
level features with the help of pre-defined anchors. The above methods
accomplish the detection procedure through two stages, including object
proposal extraction, object recognition and localization. In [12], the feature
pyramid network is exploited to further improve the detection performance of
Fast R-CNN and Faster R-CNN with the help of multi-scale feature maps. The
other pipeline for object detection implements object localization and
identification in single stage through simultaneous bounding box regression
and object classification, such as YOLO [20] and SSD [15]. The RetinaNet [13]
is also built upon the feature pyramid network, and uses dense box predictions
during the training stage. The focal loss is proposed to cope with the class
imbalance problem. The detection task has also attracted a large amount of
research interest in medical images, such as lesion detection in CT scans [25]
and cell detection in malaria images [9]. This paper targets at detecting
diseases in chest X-ray images. Practically, we propose a TBX-Att dataset to
exploit attribute information to enhance TB feature representations of disease
proposals.

###  2.2 Disease Diagnosis in Chest X-ray Images

Accurately recognizing and localizing diseases in chest X-Ray images is very
challenging because of low textural contrast, large anatomic variation across
patients, and organ overlapping. Previous works in this field mainly focus on
disease classification [17, 24, 1, 27]. Recently, the authors in [4] propose
to transfer deep models pretrained on the ImageNet dataset [5] for recognizing
pneumonia in chest X-ray images. In [22], the artificial ecosystem-based
optimization algorithm is used to select the most relevant features for
tuberculosis recognition. Based on the category activation map [29] which can
be estimated with a disease recognition network, researchers attempt to
localize disease in a weakly supervised manner [24]. In [26], the triplet loss
is used to facilitate the training of the disease classification model, and
better performance is observed in class activation maps (CAM) estimated by the
trained model. In [11, 14], multiple instance learning is employed to solve
the disease localization problem. In [28], a novel weakly supervised disease
detection model is devised on the basis of the DenseNet [8]. Two pooling
layers including a class-wise pooling layer and a spatial pooling layer are
used to transform 2-dimensional class attention maps into the final prediction
scores. The performance of these methods is still far from practical usage in
automatic diagnosis systems. In [16], a new benchmark is proposed to identify
and search potential TB diseases, however, it lacks attribute information that
has been proven to pose a great effect on the practice diagnosis of medical
experts.

Figure 1:  Illustration of the proposed Attribute Relational Reasoning
Network.

##  3 Methods

###  3.1 Overview

The method proposed in this paper mainly addresses the weakly supervised
tuberculosis detection from chest X-rays. The pipeline of our proposed method
is highlighted in Fig. 1, which include three parts: (1) the feature pyramid
network is utilized for extracting multi-scale tuberculosis feature maps; (2)
the attribute classifier, which employs the attribute supervised information
to generate multi-scale attribute feature maps; (3) the feature interaction
module takes the attribute feature map to establish the attribute prompt to
obtain a more representative feature representation of tuberculosis. Compared
with the original Feature Pyramid Network used for object detection, another
branch is added to generate the corresponding multi-scale features for
attribute classification, which can prompt a guidance to extract
representative features for detection.

###  3.2 Attribute Feature Representation

Figure 2:  Attribute feature representation. (a) shows the process of
attribute feature extraction. (b) shows the fusion of multi-scale attribute
features and the classification.

The detection branch is designed to utilize multi-scale features to generate
foreground detection result. The backbone network used for abstracting feature
representation always generate multi-scale features. Specifically, for ResNets
we use the feature activation output by each stage’s last residual block. We
denote the output of these blocks as $\\{C_{2},C_{3},C_{4},C_{5}\\}$, which
have strides of $\\{4,8,16,32\\}$ pixels with respect to the input image.
Correspondingly, multi-scale attribute features as shown in Fig. 2 are
constructed in our method for both the task of classification and feature
interaction described in the next section.

Given the input $C_{i}$ from the i-th stage’s last residual block, we simply
attach a $1\times 1$ convolution layer to produce a feature map having
$C_{a}\times N_{a}$ channels, where $N_{a}$ represents the number of kinds of
attributes and $C_{a}$ represents the channel dimension of each attribute
feature map. Then two group convolution layers are used to generate the
distinguishable feature maps for each attribute individually. For each group
convolution layer, we set the number of groups as $N_{a}$ and expect each
group to represent a kind of attribute. For better taking advantage of
information flowing across attributes’ feature maps, channel shuffle operation
is added in the last group convolution layer. The specially designed Attribute
Block used for extracting attribute features is shown in Fig. 2(a).

After obtaining the features of each scale($A_{2}$ to $A_{5}$ in Fig. 2(b),
the Global Average Pooling operation and feature fusion by concatenation are
used to get the final feature vectors used for attribute classification.

###  3.3 Feature Interaction

Multi-head Cross-Attention Given two feature maps X, $\textbf{Y}\in R^{C\times
H\times W}$, where $H$,$W$ are the spatial height, width and $C$ is the number
of channels. One $1\times 1$ convolution is used to project X to query
embedding $\textbf{Q}_{i}$, another two $1\times 1$ convolutions are adopt
similarly to project Y to key embedding $\textbf{K}_{i}$ and value embedding
$\textbf{V}_{i}$, where $i$ represents the i-th head. The dimension of
embedding in each head is $d$, and the number of heads is $N=C/d$. The
$\textbf{Q}_{i}$, $\textbf{K}_{i}$, $\textbf{V}_{i}$ are then flattened and
transposed into sequences with the size of $n\times d$, where $n=H\times W$.
The output of the cross-attention (CA) in i-th head is a scaled dot-product:

$CA_{i}(\textbf{Q}_{i},\textbf{K}_{i},\textbf{V}_{i})=SoftMax(\frac{\textbf{Q}_{i}\textbf{K}_{i}^{T}}{\sqrt{d}})\textbf{V}_{i}.$
(1)

To reduce the large computational complexity of CA on high resolution feature
maps, $\textbf{K}_{i}$ and $\textbf{V}_{i}$ are down-sampling by the ratio of
$s$, resulting smaller feature maps with the size of
$d\times\frac{H}{s}\times\frac{W}{s}$. In our implementation, $s$ is set to 16
in the first stage, and gradually reduce by 2 times. Finally, in a $N$-head
attention situation, the output after multi-head cross-attention (MCA) is
calculated as follows:

$MCA(\textbf{X},\textbf{Y})=\Phi(Concat(CA_{1},CA_{2},...,CA_{N}))$ (2)

where $\Phi(\cdot)$ is a liner projection layer that weights and aggregates
the feature representation of all attention heads.

####  3.3.1 Attr-Attr Attention ($A^{2}$-Attn)

Figure 3:  Cross-attention. (a) shows the multi-head cross attention among
different kinds of attribute features. (b) shows cross attention between
attribute features and TB features.

Given the input
$\textbf{X}=\\{\textbf{F}_{1},\textbf{F}_{2},...,\textbf{F}_{N_{a}}\\}$ from
the i-th scale of attribute feature maps, attention modules are used to
reproduce more representative features by implicitly mining the relationship
among $N_{a}$ kinds of attributes. To build the attention module, we firstly
concatenate $\textbf{F}_{i}$ channel-wisely, and then linearly projected to
get Y, which has the same channel dimension of $\textbf{F}_{i}$. As shown in
Fig. 3(a), the attention module takes $\\{\textbf{F}_{i},i=1,...,N_{a}\\}$ as
queries and uses Y to generate key and value embeddings, so the finally
aggregated feature $\textbf{F}_{i}^{{}^{\prime}}$ can be obtained by
$\textbf{F}_{i}^{{}^{\prime}}=MCA(\textbf{F}_{i},\textbf{Y}$). The operation
of $A^{2}$-Attn can be applied to attribute features of each scale ($A_{2}$ to
$A_{5}$ in Fig. 2(b)) to get more representative features.

####  3.3.2 Attr-TB Attention (AT-Attn)

Given attribute features $\\{\textbf{F}_{i},i=1,...,N_{a}\\}$ and TB feature
$\textbf{F}_{tb}$, we expect to design a attribute and TB attention module
(AT-Attn) to enhance the TB representation with attribute prompt. We firstly
reaggregate attribute features to get X by taking advantage of the similarity
between attribute features and TB feature. Specifically, GAP operation
followed by a linear projection is performed to get attribute vectors
$\\{\textbf{A}_{i},i=1,...,N_{a}\\}$ and TB vector B. The similarity score
$s_{i}$ is calculated by $\textbf{A}_{i}\cdot\textbf{B}$, and then X can be
obtained by $\textbf{X}=\sum s_{i}\textbf{F}_{i}$. As shown in in Fig. 3(b),
the attention module takes X as query and use $\textbf{F}_{tb}$ to generate
key and value embeddings. Finally, the output of the AT-attn is obtained as
follows:

$\textbf{F}_{tb}^{{}^{\prime}}=\textbf{F}_{tb}+Norm(MCA(\sum
s_{i}\textbf{F}_{i},\textbf{F}_{tb}))$ (3)

###  3.4 Loss Function

During the training, the classification branch and the detection branch are
jointly trained, which can be seen as a type of multi-tasking learning.
$Loss_{det}$ is used to optimize the detection branch, and $Loss_{cls}$ is
performed to train the classification branch. For $Loss_{cls}$, we apply the
most commonly used binary cross-entropy loss. For $Loss_{det}$, the same loss
function is used as in [21], which contains an anchor category prediction item
and a bounding box regression item. The final loss function is the sum of
$Loss_{det}$ and $Loss_{cls}$, and $\lambda$ is used for balancing the
influence of two loss terms.

$Loss=Loss_{det}+\lambda Loss_{cls}$ (4)

###  3.5 Experimental Analysis

####  3.5.1 Dataset

The collection of attribute-based TB X-rays faces many difficulties [16]. The
attribute information is so hard to obtain due to its unclear pathological
structure. To reduce these problems, we cooperate with one of the most
important hospitals in China to collect X-rays. We also combine the TBX11K
dataset [16] to construct a totally new attribute-based TB dataset. The TBX11K
dataset consists of 11200 X-rays, including 5000 healthy cases, 5000 sick but
non-TB cases, and 1200 cases with manifestations of TB. The dataset includes 8
types of sign: Fibrotic Streaks, Pulmonary Consolidation, Diffuse Nodules,
Pulmonary Cavitation, Atelectasis, Multiple Nodules, Pleural Effusion, and
Pulmonary Tuberculosis. Note that all X-ray images are resized and have the
resolution of about 512 $\times$ 512\. Two senior radiologists participated in
the annotation independently following a list of attributes. If an attribute
is present, it is annotated as 1 and otherwise 0, resulting in a sparse vector
of attributes. The sample was retained only if the two radiologists labeled
the same. The attribute dataset consists of a toal of 2000 X-ray images,
including a training set of 1700 images and a validation set of 300 images.
For the combination with the existing TBX11K dataset, annotations from the
attribute dataset are used to train the attribute classification branch, while
the TBX11K dataset consists of TB localization information which provides the
supervision for training the detection branch.

####  3.5.2 Settings

We adopt Faster RCNN as baseline, in addition, a single scale feature
interaction is firstly used to explore the influence of cross attribute
attention. Then a multi-scale feature interaction way is adopt to further
enhance the attention function. For all experiments, each model is
independently trained five times with randomly initialized weights and
correspondingly 5 validations are performed for each model. The overall
performance of a model is assessed with several commonly used metrics,
including the mean and standard deviation of accuracy, F-score and mAP. All
models are trained for 60 epoches from scratch using PyTorch [18] on NVIDIA
Titan X pascal GPUs while Adam [10] being the optimizer with the initial
learning rate set to 1e-3, which is reduced by a factor of 10 after every 20
epochs. The weight decay is set to 1e-4. We adopt ResNet-50 as the backbone
model. The batch size is 8 on a single GPU. Different from [16], we train the
two branches simultaneously rather than dividing it into two stages.

Table 1: Performance Comparison of Tuberculosis Diagnosis Models on the TB-
Xatt Dataset.

Methods | Detector | Major Component | Results(%) (mean$\pm$standard deviation)
---|---|---|---
GroupConv |  $A^{2}$-Attn | AT-Attn | F-score | Accuracy | mAP
Baseline | Two-stage Model | $\boxtimes$ | $\boxtimes$ | $\boxtimes$ | 29.24$\pm$0.76 | 88.08$\pm$0.12 | 17.10$\pm$0.11
SingleScale | Two-stage Model | $\boxtimes$ | $\checkmark$ | $\checkmark$ | 33.70$\pm$0.63 | 92.08$\pm$0.71 | 18.20$\pm$0.17
$\checkmark$ | $\boxtimes$ | $\checkmark$ | 32.70$\pm$0.72 | 91.08$\pm$0.12 | 17.20$\pm$0.16
$\checkmark$ | $\checkmark$ | $\boxtimes$ | 32.42$\pm$0.12 | 91.03$\pm$0.64 | 17.09$\pm$0.13
|  | $\checkmark$ | $\checkmark$ | $\checkmark$ | 34.24$\pm$0.63 | 93.15$\pm$0.21 | 19.00$\pm$0.12
MultiScale | Two-stage Model | $\boxtimes$ | $\checkmark$ | $\checkmark$ | 33.13$\pm$0.14 | 94.01$\pm$0.21 | 18.79$\pm$0.32
$\checkmark$ | $\boxtimes$ | $\checkmark$ | 34.12$\pm$0.23 | 91.82$\pm$0.33 | 17.92$\pm$0.42
$\checkmark$ | $\checkmark$ | $\boxtimes$ | 34.22$\pm$0.11 | 92.01$\pm$0.24 | 18.01$\pm$0.22
|  | $\checkmark$ | $\checkmark$ | $\checkmark$ | 39.37$\pm$0.12 | 94.61$\pm$0.11 | 19.20$\pm$0.10

####  3.5.3 Results

On our new proposed TBX-Att dataset for weakly tuberculosis detection, we
compared our attribute relational reasoning with a two-stage faster RCNN
model. We focus on the impacts of three carefully designed components on both
attribute classification and TB detection in experiments. As shown in Table 1,
when the three components are all used, the highest performance is obtained
and exceeds baseline by a large margin. Our proposed method under multi-scale
setting achieves an accuracy of 94.61%, a F-score of 39.37% and a mAP of
19.20%, and outperforms the baseline by 10.13%, 6.53% and 2.1% respectively.
Any pairwise combination will reduce the performance, which can demonstrate
the effectiveness of each component implicitly. Some interesting phenomenon
that the feature interaction can gain improvements on both tasks (not just TB
detection) is found in Table 1. Moreover, the results verify that the two
tasks are complementary, where the higher accuracy of attribute classification
leads to higher TB detection performance. And the above findings also confirm
our original intention, which takes advantage of TB attribute information to
guide the effective extraction of features for TB detection.

##  4 Conclusion

In this paper, we present a TBX-Att dataset with attribute labels to expand
the existing TBX11K dataset for weakly tuberculosis detection, considering
radiologists usually use clinical features like attribute information for
diagnosis. Moreover, a multi-scale feature interaction model is devised to
enhance TB feature representations under the guidance of relational knowledge
reasoning. We find out that with the attribute information assistance, TB
classification and detection are more easily to achieve a better performance.
We hope this dataset can further inspire not only researchers but also medical
experts.

####  4.0.1 Acknowledgements

This work is supported by the National Natural Science Foundation of China
(Grant Nos. 62141605, 62106248, U21B2048).

## References

  * [1] Aviles-Rivero, A.I., Papadakis, N., Li, R., Sellars, P., Fan, Q., Tan, R.T., Schönlieb, C.B.: Graphxnet chest x-ray classification under extreme minimal supervision. arXiv:1907.10085 (2019)
  * [2] Brestel, C., Shadmi, R., Tamir, I., Cohen-Sfaty, M., Elnekave, E.: Radbot-cxr: Classification of four clinical finding categories in chest x-ray using deep learning. In: Proc. Med. Imag. Deep. Learn. pp. 1–9 (2018)
  * [3] Cai, J., Lu, L., Harrison, A.P., Shi, X., Chen, P., Yang, L.: Iterative attention mining for weakly supervised thoracic disease pattern localization in chest x-rays. In: Frangi, A.F., Schnabel, J.A., Davatzikos, C., Alberola-López, C., Fichtinger, G. (eds.) Medical Image Computing and Computer Assisted Intervention – MICCAI 2018. vol. 11071, pp. 589–598. Springer, Cham (2018). https://doi.org/10.1007/978-3-030-00934-2_66
  * [4] Chouhan, V., Singh, S.K., Khamparia, A., Gupta, D., Tiwari, P., Moreira, C., Damaševičius, R., De Albuquerque, V.H.C.: A novel transfer learning based approach for pneumonia detection in chest x-ray images. Applied Sciences 10(2), 559 (2020)
  * [5] Deng, J., Dong, W., Socher, R., Li, L.J., Li, K., Fei-Fei, L.: Imagenet: A large-scale hierarchical image database. In: 2009 IEEE conference on computer vision and pattern recognition. pp. 248–255. Ieee (2009)
  * [6] Girshick, R.: Fast r-cnn. In: Proc. IEEE Int. Conf. on Comput. Vis. (ICCV). pp. 1440–1448 (2015)
  * [7] Girshick, R., Donahue, J., Darrell, T., Malik, J.: Rich feature hierarchies for accurate object detection and semantic segmentation. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 580–587 (2014)
  * [8] Huang, G., Liu, Z., Van Der Maaten, L., Weinberger, K.Q.: Densely connected convolutional networks. In: Proc. IEEE Conf. Comput. Vis. Pattern Recognit. (CVPR). pp. 4700–4708 (2017)
  * [9] Hung, J., Carpenter, A.: Applying faster r-cnn for object detection on malaria images. In: Proc. IEEE Conf. Comput. Vis. Pattern Recognit. (CVPR) Wksp. pp. 56–61 (2017)
  * [10] Kingma, D.P., Ba, J.: Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980 (2014)
  * [11] Li, Z., Wang, C., Han, M., Xue, Y., Wei, W., Li, L.J., Fei-Fei, L.: Thoracic disease identification and localization with limited supervision. In: Proc. IEEE Conf. Comput. Vis. Pattern Recognit. (CVPR). pp. 8290–8299 (2018)
  * [12] Lin, T.Y., Dollár, P., Girshick, R., He, K., Hariharan, B., Belongie, S.: Feature pyramid networks for object detection. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 2117–2125 (2017)
  * [13] Lin, T.Y., Goyal, P., Girshick, R., He, K., Dollár, P.: Focal loss for dense object detection. In: Proc. IEEE Int. Conf. on Comput. Vis. (ICCV). pp. 2980–2988 (2017)
  * [14] Liu, J., Zhao, G., Fei, Y., Zhang, M., Wang, Y., Yu, Y.: Align, attend and locate: Chest x-ray diagnosis via contrast induced attention network with limited supervision. In: Proc. IEEE Int. Conf. on Comput. Vis. (ICCV). pp. 10632–10641 (2019)
  * [15] Liu, W., Anguelov, D., Erhan, D., Szegedy, C., Reed, S., Fu, C.Y., Berg, A.C.: Ssd: Single shot multibox detector. In: Leibe, B., Matas, J., Sebe, N., Welling, M. (eds.) Computer Vision – ECCV 2016. pp. 21–37. Springer, Cham (2016)
  * [16] Liu, Y., Wu, Y.H., Ban, Y., Wang, H., Cheng, M.M.: Rethinking computer-aided tuberculosis diagnosis. In: Proceedings of the IEEE/CVF conference on computer vision and pattern recognition. pp. 2646–2655 (2020)
  * [17] Mohd Noor, N., Mohd Rijal, O., Yunus, A., Mahayiddin, A.A., Gan, C.P., Ong, E.L., Abu Bakar, S.A.R.: Texture-based statistical detection and discrimination of some respiratory diseases using chest radiograph. In: Advances in Medical Diagnostic Technology, pp. 75–97. Springer, Singapore (2014)
  * [18] Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., et al.: Pytorch: An imperative style, high-performance deep learning library. In: Advances in neural information processing systems. pp. 8026–8037 (2019)
  * [19] Qin, C., Yao, D., Shi, Y., Song, Z.: Computer-aided detection in chest radiography based on artificial intelligence: a survey. Biomed. Engg. OnLine 17(1), 113 (Aug 2018)
  * [20] Redmon, J., Divvala, S., Girshick, R., Farhadi, A.: You only look once: Unified, real-time object detection. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 779–788 (2016)
  * [21] Ren, S., He, K., Girshick, R., Sun, J.: Faster r-cnn: Towards real-time object detection with region proposal networks. In: Advances in neural information processing systems. pp. 91–99 (2015)
  * [22] Sahlol, A.T., Abd Elaziz, M., Tariq Jamal, A., Damaševičius, R., Farouk Hassan, O.: A novel method for detection of tuberculosis in chest radiographs using artificial ecosystem-based optimisation of deep neural network features. Symmetry 12(7), 1146 (2020)
  * [23] Uijlings, J.R., Van De Sande, K.E., Gevers, T., Smeulders, A.W.: Selective search for object recognition. Int. J. Comput. Vis. 104(2), 154–171 (2013)
  * [24] Wang, X., Peng, Y., Lu, L., Lu, Z., Bagheri, M., Summers, R.M.: Chestx-ray8: Hospital-scale chest x-ray database and benchmarks on weakly-supervised classification and localization of common thorax diseases. In: Proc. IEEE Conf. Comput. Vis. Pattern Recognit. (CVPR). pp. 2097–2106 (2017)
  * [25] Yan, K., Wang, X., Lu, L., Summers, R.M.: Deeplesion: automated mining of large-scale lesion annotations and universal lesion detection with deep learning. J. Med. Imag. 5(3), 036501 (2018)
  * [26] Zhang, C., Chen, F., Chen, Y.Y.: Thoracic disease identification and localization using distance learning and region verification. arXiv:2006.04203 (2020)
  * [27] Zhao, G., Feng, Q., Chen, C., Zhou, Z., Yu, Y.: Diagnose like a radiologist: Hybrid neuro-probabilistic reasoning for attribute-based medical image diagnosis. IEEE Transactions on Pattern Analysis and Machine Intelligence pp. 1–1 (2021)
  * [28] Zhou, B., Li, Y., Wang, J.: A weakly supervised adaptive densenet for classifying thoracic diseases and identifying abnormalities. arXiv:1807.01257 (2018)
  * [29] Zhou, B., Khosla, A., Lapedriza, A., Oliva, A., Torralba, A.: Learning deep features for discriminative localization. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 2921–2929 (2016)