Evaluating online filtering algorithms to enhance dynamic multispectral optoacoustic tomography.

Multispectral optoacoustic tomography (MSOT) is an emerging imaging modality, which is able to capture data at high spatiotemporal resolution using rapid tuning of the excitation laser wavelength. However, owing to the necessity of imaging one wavelength at a time to the exclusion of others, forming a complete multispectral image requires multiple excitations over time, which may introduce aliasing due to underlying spectral dynamics or noise in the data. In order to mitigate this limitation, we have applied kinematic α and α β filters to multispectral time series, providing an estimate of the underlying multispectral image at every point in time throughout data acquisition. We demonstrate the efficacy of these methods in suppressing the inter-frame noise present in dynamic multispectral image time courses using a multispectral Shepp-Logan phantom and mice bearing distinct renal cell carcinoma tumors. The gains in signal to noise ratio provided by these filters enable higher-fidelity downstream analysis such as spectral unmixing and improved hypothesis testing in quantifying the onset of signal changes during an oxygen gas challenge.

Introduction

Imaging has become a powerful tool in diverse fields such as medical diagnostics, high-throughput screening, and remote sensing [[1], [2], [3], [4], [5], [6], [7], [8], [9], [10], [11]]. The imaging toolkit is becoming progressively more sophisticated, moving beyond spatial discrimination with additional dimensions: Spectral imaging allows the capture of multiple imaging channels, in the hope of using the spectral fingerprints of underlying components to learn more about the system under study [1,[7], [8], [9], [10],[12], [13], [14], [15], [16], [17], [18], [19]]; and dynamic imaging captures many single images over time in order to resolve changes in the observed system [12,13,[20], [21], [22], [23], [24]]. Imaging methods are subject to inherent trade-offs between spatial, spectral, and temporal resolution, as well as the discrimination of signal versus noise. The quality of data is fixed once the data have been acquired by instrument hardware, but it is possible to take advantage of the inherent coherence and structure within a dataset in order to mitigate the effects of noise and finite sampling.

Multispectral optoacoustic tomography (MSOT) is an imaging modality, characterized by the use of optical excitation and ultrasonic detection, which is able to capture spectral images at high spatial and temporal resolution [7,17,20,22]. MSOT intrinsically acquires data based on sequential wavelength excitation such that images are captured one channel at a time. There may therefore be substantial spectral distortion due to aliasing in imaging scenarios where the rate of spectral dynamics is comparable to the multispectral sampling rate. Further exacerbating aliasing is the presence of noise, which may vary spatially, temporally, and spectrally, resulting in the reduction of measurement fidelity.

One experimental motif which exemplifies these difficulties is the gas-breathing challenge, wherein the inhaled gas of a subject under observation is varied in order to induce a change in the subject’s blood oxygen saturation (SO2) [25,26]. In order to determine or estimate the oxygen saturation by MSOT, images must be acquired at multiple wavelengths and spectrally unmixed to provide oxy- and deoxy-hemoglobin images, whereupon the ratio of oxy-hemoglobin to total hemoglobin provides an estimate of the oxygen saturation of blood. The propagation of image noise through this process can result in extremely noisy SO2 estimates, limiting the amount of useful information which can be extracted from noisy data.

It is therefore desirable to devise a scheme to mitigate the corrupting effects of noise while still accurately capturing dynamic phenomena. Similarly, it is desirable to be able to perform such filtering in an online manner, as data are acquired, such that the filtering may be implemented for real-time imaging or efficient post-acquisition analysis. This last criterion imposes a strong restriction on the complexity of any algorithm, since many photoacoustic imaging systems gather data at rates exceeding 10 Hz (.

We briefly establish general notation for equations throughout this work: is the current time index, refers to the states of interest, here the pixel values for the spectral image; is the acquired data at time , is the estimate of at time k considering all data up to and including the data at time , is the so-called prior residual, and is the estimate of after integration of the data acquired at time point .

A simple, efficient approach to recursively estimating the state of a multispectral image is to perform online averaging on a per-channel basis according to the following update equations:

For stationary distributions, online averaging performs exactly as the usual average, converging to the true mean as , i.e., the variance of sampling noise is proportional to the square root of the number of samples averaged [27]. For non-stationary distributions such as those encountered in dynamic imaging, will not, in general, converge to the true mean even after acquisition of an infinite number of samples. In the context of a gas challenge, will be heavily biased towards the initially acquired data and will not effectively follow dynamics.

To accommodate nonstationary distributions, the term may be replaced with a constant . This gives the averaging process a time-invariant form which may be analyzed as a recursive filter, often referred to as the exponentially weighted moving average or filter:

Observation indicates that varying between 0 and 1 biases the estimate of towards the previous estimate and the new data, respectively. converges to the same structure as the online average as the sampling time tends to infinity. Conversely, ignores all previous data in favor of the new data, but without any noise rejection. This trade-off between noise rejection and dynamical sensitivity is irrelevant if there are no underlying dynamics, but for values of the estimate will always lag the true mean of dynamic systems [28]. The limiting case of results in the filter taking on the form of a sliding-window (SW) filter, where old data is entirely overwritten by new data as it is acquired:

Clearly, a simple filter is unable to fully accommodate nonstationary data and is therefore generally inappropriate for dynamic imaging. Many more sophisticated algorithms such as the broad class of Bayesian filters can incorporate large amounts of prior knowledge, including the presence of dynamics, and have exceptional performance in tracking the state of dynamical systems [28,29]. However, Bayesian filters can require large amounts of data to undergo computationally intensive processing and are particularly challenging to apply to spectral images due to the large size of associated datasets.

Though full Bayesian filters are computationally prohibitive, they may be simplified by assuming Gaussian noise statistics and linear dynamics. These so-called Kalman filters operate in a two-stage prediction-correction cycle, wherein the prediction phase propagates the state estimate in time according to the dynamical model and the correction phase incorporates new data observations into the state estimate. Kalman filters have been shown to be useful in integrating multi-domain measurements, including multi-channel images, due to their inherently dynamical nature [[30], [31], [32], [33], [34]]. Nevertheless, the size of the covariance matrices necessary for the operation of Kalman filters requires offline processing or substantial computing resources.

If the assumption is made that the Kalman filter is operating at steady state, it degenerates to a form which is temporally invariant. Further assumption that the state of interest has a known second-order process variance (i.e., that the process noise is in the acceleration term of a Newtonian kinematic model) yields a family of recursive filters indexed by the order of their Taylor series expansion, with associated constant coefficients [35,36]. Here, we describe the form of the 1st and 2nd order filters, the former of which is the previously-described filter and the latter correspondingly named the filter.

As before, the state of interest is denoted by , with first temporal derivative given by . The measurement noise variance, describing the ‘noisiness’ of each measurement, is denoted by with units determined by the measurement device (here arbitrary units a.u.), and the process noise variance, describing the volatility of the underlying dynamics, as with units (a.u./s2). We assume that each pixel evolves independently and thus the filter is applied to each pixel of the acquired image separately.

The filter coefficients may be optimally found using prior knowledge of , , and via a single scalar parameter, the tracking index:

We assume that the values of , , and are homogenous across the entire image such that each image has a single tracking index. In practice, these parameters may vary between channels and over time, which we address in Section 4.

The filtering equations for an filter are identical to the exponentially weighted moving average described in Eqs. 1 and 3. As the filter tracks only position information, the prediction step is degenerate and only the correction phase updates the state estimate:

Prediction

Correction

The filtering equations for an filter propagate the state estimate over time by assuming constant-velocity dynamics and are:

Prediction

Correction

The tracking index formalism enables analytic determination of the and coefficients (Appendix A) in an optimal manner, and requires knowledge of only three parameters , , and , which may themselves be determined from prior knowledge of the system, or updated recursively by examining the prior and posterior residuals of the filtering process [37]. In environments where the sampling rate is inhomogeneous, such as when oversampling each channel before sampling another, is variable. This requires that the tracking index and thus the and coefficients must be updated at each time point.

We compare the and filters against the sliding window (SW) filter equivalent to : A circular buffer of data is maintained and updated on a per-channel per-time-point basis as data are acquired, following a ‘sample and hold’ pattern. There is thus no incorporation of previous data into the current estimate, nor is there any accommodation for dynamic evolution of the state.

Methods

Numerical phantom

We developed a dynamic, multispectral numerical phantom to act as a known ground-truth for testing the ability of the various filters to track spectral changes while rejecting noise. The full process of creating the phantom is described in Appendix B and illustrated in Fig. 1. Briefly, a Shepp-Logan phantom [38] was modified to include time-varying concentrations of endmembers with known spectra. Sampling the phantom image at specific wavelengths at each point in time provided a sample train to emulate the acquisition of one channel at a time. We directly sampled the phantom, rather than simulating the photoacoustic imaging process, to limit the possibility of introducing implementation and modeling errors. After sampling, zero-mean white Gaussian noise with standard deviation was added to each image, corresponding to an SNR of approximately 20. The dynamics of the phantom were designed to mimic a gas challenge due to the variety of states and rates of change which can be present in a single time course.

Fig. 1

Workflow for creating numerical phantom. A) A Shepp-Logan phantom was created with ellipses defined by various constant values for total hemoglobin [Hbtot] and time-varying SO2. The time-varying images of [Hb] and [HbO2] were then calculated (B) and their spectra used to create time-varying spectral images (C). The spectral image was then sampled at predetermined wavelengths to create a time course of single-wavelength images (D), which were used as input for the filtering methods described in this work.

Experimental data

Three mice from ongoing studies were used for in vivo validation, with all data acquired using an MSOT InVision 256-TF photoacoustic imaging system [39] (iThera Medical GmbH, Munich, Germany). The water bath used for imaging was maintained at a temperature of 35 C for all experiments. All animal work was conducted in accordance with the UT Southwestern Institutional Animal Care and Use Committee guidelines as well as all superseding federal guidelines. Two distinct wavelength sampling patterns were used to demonstrate the applicability of the filters to general acquisition conditions.

The first mouse was a female BALB/c mouse (UT Southwestern Animal Breeding Core) with a RENCA-luc [40] mouse kidney tumor implanted in the right kidney following published protocols [41], with data acquired 28 days after implantation. Prior to imaging, the mouse was shaved and any remaining hair was removed from the animal using depilation cream. Five wavelengths were sampled [λ = 700, 730, 756, 800, 850 nm] to correspond with major features of the hemoglobin spectrum, including local spectral peaks and the isosbestic point. Each wavelength was oversampled 5 times and averaged prior to image reconstruction, leading to an overall acquisition time of 2.5 s per multispectral sampling cycle. All images were acquired at the same axial cross-section. During imaging, the anesthetized mouse underwent a gas challenge from air-oxygen-air, with 10.5 min on air, 10 min on oxygen, followed by another 5 min on air, while breathing 2.0 % isoflurane at a rate of 2 L/min.

The remaining two mice were female Hsd:Athymic Nude-Foxn1nu mice (Envigo) with XP373 human renal cell carcinoma tumor xenografts [42] implanted in the right kidney following published protocols [41], with two imaging sessions acquired for Mouse #1 at 31 and 36 days, and one session at for Mouse #2 at 31 days after tumor implantation. Prior to imaging, hair was removed from the animal using depilation cream. Seven wavelengths were sampled sequentially [λ = 680, 734, 757, 798, 850, 922, 980 nm] to correspond with major features of the hemoglobin spectrum, including local spectral peaks and the isosbestic point. Each wavelength was sequentially oversampled 3 times and averaged prior to image reconstruction, leading to an overall acquisition time of 2.1 s per multispectral sampling cycle. All images were acquired at the same axial cross-section. During imaging, each anesthetized mouse underwent a gas breathing challenge from air-oxygen-air, with 5 min on each gas, while breathing 2.0 % isoflurane at a rate of 2 L/min.

Image reconstruction

After acquisition, each single-wavelength image was preprocessed by correcting for laser power variation and deconvolving the signal impulse response using Wiener deconvolution with a regularization constant of 0.001. Each image was then reconstructed using a curve-driven model matrix as a forward model [43] and a non-negative accelerated projected conjugate gradient method as the inverse solver [44].

Multispectral state filtering

All analysis was performed using code developed in MATLAB 2017b, with production code available at https://git.biohpc.swmed.edu/PIRL/kalata-filters-matlab. Archived version is available at https://doi.org/10.5281/zenodo.3836682.

The single-wavelength images from the numerical and experimental datasets were used as input to each of the SW, , and filters. For each method, an multispectral image with pixels per channel and channels was maintained, to store each filter’s estimate of the multispectral image. At each point in time, all channels of the entire multispectral image were updated according to the prediction step for each method, using either Eq. 6 for the SW and filters, or Eqs. 8 and 9 for the filter. Following the prediction step, a single-channel image, acquired at wavelength , was loaded into memory, corresponding to in Eq. 1. This image was used to update channel in each filter’s multispectral state estimate, using either Eq. 7 for the SW and filters or Eqs. 10 and 11 for the filter.

The numerical phantom data were filtered using , , the RENCA-luc data with , , and the XP373 data with ,

Post-processing

At each time point, the multispectral images estimated by each filtering method were spectrally unmixed using an unconstrained linear model to yield images of oxy- and deoxy-hemoglobin. The total hemoglobin and oxygen saturation in each pixel were estimated as:

We note that the resultant values are estimates, and are thus indicated by the superscript MSOT as used by others [26]; truly quantitative measures of endmember concentration would require additional processing using, e.g., fluence correction [18,19,45,46] or Bayesian methods [47]. In general, the filtering methods presented here cannot correct for errors in quantitation due to modelling errors.

Quantification methods

Image quality

The deviation of each estimated multispectral image from the known ground truth image was calculated using two metrics: the mean-squared error (MSE) and the structural similarity index (SSIM).

The MSE was used to examine the effects of each approach on noise rejection, but includes some deviation due to differences in the image mean. It is calculated as:

The SSIM was used to examine images on a holistic basis, accounting for structure, noise, and constant value offsets. It is calculated as:where and represent the mean of each image, and represent the variance of each image, and represents the covariance between the two images. and are positive constants chosen to improve numerical stability and are related to the dynamic range of the images.

ROI analysis

For the phantom data, a 5 × 5 region of interest (ROI) was defined in the same location for all images within ellipse (e) (Appendix B, Table B1), and the average value within that region determined for all time points for each channel. This was performed for all filtering methods, for each channel of the multispectral image, as well as for each of the [Hb]MSOT, [HbO2]MSOT, [Hbtot]MSOT, and SO2 MSOT images, to yield time courses of each parameter.

Table B1

Parameters for dynamic multispectral phantom.

Ellipse	Center (x,y)	Major Axis	Minor Axis	Theta	Hbtot	c0	c1	c2	kresponse
a	(0,0)	0.69	0.92	0	1	0.4	0.6	0.4	0.2
b	(0, −0.0184)	0.6624	0.874	0	0.8	1.0	0.5	0.0	0.3
c	(0.22,0)	0.11	0.31	−18°	0.2	0.2	0.5	0.1	0.02
d	(−0.22, 0)	0.16	0.41	18°	0.2	0.5	0.5	0.6	0.04
e	(0, 0.35)	0.21	0.25	0	0.1	0.2	0.9	0.2	0.2
f	(0, 0.1)	0.046	0.046	0	0.1	0.1	0.1	0.9	0.1
g	(0, −0.1)	0.046	0.046	0	0.1	0.9	0.8	0.7	0.05
h	(−0.08, −0.605)	0.046	0.023	0	0.1	1.0	1.0	1.0	0.4
i	(0, −0.605)	0.023	0.023	0	0.1	0.7	0.6	0.2	0.2
j	(0.06, −0.605)	0.023	0.046	0	0.1	0.9	0.2	0.9	0.1

Manually selected ROIs were defined over the tumor and spine for each mouse and used as the ROIs for both single-channel data and the unmixed parameters. Within each ROI, for each filter, for each parameter, at each time point, a histogram was calculated using 256 equally-spaced bins, with lower and upper limits defined by the 0.1st and 99.9th percentile of the values for that parameter in that ROI for all time points. Each one-dimensional histogram was then used as a column of pixels in a dynamic histogram image. Three time points for each tumor type were chosen prior to, during, and after the transition period from air to oxygen. For each condition, the histograms from 11 time points centered on those selected time points were averaged together and used as representative time slices of the dynamic histogram image to improve visualization.

Further analysis was performed on each filter to compare the statistical effects of each, and this is described in Appendix C.

Results

Phantom data

When applied to the simulated data, both and filters effectively suppressed noise in the multispectral images as compared to the SW filter. As seen in Fig. 2, the filter rejected a large amount of noise from the signal, but failed to follow the rapid dynamics of the simulated gas challenge, causing visible artifacts (Fig. 2A, bottom right). By comparison, the filter both rejected noise and accurately tracked dynamic changes throughout the image. The relative performance of each of the filtering methods can be seen in the single-channel ROI time course and quality metric time courses. Overall, the filter performed extremely well unless challenged by sudden dynamic changes. Due to the filter’s intrinsic time lag, there was a substantial resettling time before it returned to pre-challenge performance. The properties of each filter are retained upon spectral unmixing of the multispectral images, and still further when calculating the derived [Hbtot]MSOT and SO2MSOT images as seen in Fig. 3. The filter demonstrated superior performance on the unchanging [Hbtot]MSOT channel, but once again suffered from a time lag following dynamic changes. These observations are validated by the time courses of error metrics in Fig. 4. Notably, the unmixed images using the SW were extremely noisy, with the [HbO2]MSOT SSIM dropping below 0.8 towards the end of the challenge. The filter is able to achieve excellent noise rejection at baseline, yielding small MSE and SSIM very close to 1, but during transitions its performance fell below that of the SW.

Fig. 2

Effect of each filter on dynamic images in numerical phantom. A) Single-wavelength images (left) and their differences (right) from the known ground truth image (left, top) at the final time point of the simulated gas challenge. A clear reduction in noise is seen for both the and filters when compared to the SW filter. The filter suffers from a visible artifact due to its intrinsic time lag, as indicated by the yellow arrow. B) Image showing ROI for signal intensity calculation in C. C) Mean signal intensity over the yellow ROI in B for each time point during challenge. Conspicuous lag is visible in the filtered data. Red: SW, blue: , green: black: Ground truth. D) Mean squared error of entire multispectral image at each point in time. Degradation of filter performance is clearly seen accompanying gas challenge changes. E) Structural similarity of entire multispectral image at each point in time, matching D.

Fig. 3

Effect of each filter on time courses of derived parameters. A) Derived images from filtered data for each filter. The effects of noise in the multispectral images are magnified in the unmixed data. Filtering prevents noise from being propagated to this stage of analysis, leading to qualitatively clearer images for both the and filters. B) Difference of each derived parameter image from the corresponding ground truth (A, top row). Both and demonstrate decreased noise relative to SW, but the filter’s lag has resulted in an inaccurate unmixing result, overestimating [HbO2] and underestimating [Hb] as indicated by yellow arrows. C) Time courses of each parameter for different filters. Red: SW, blue: , green: black: Ground truth. Though dynamic changes cause the filter to deviate from ground truth, it has excellent performance in filtering the constant [Hbtot] channel.

Fig. 4

Image quality metrics over time for unmixed parameter images in phantom. Red: sliding window filter, green: , blue: . As in the multispectral case, the filter achieves excellent noise reduction in steady-state contexts, such as the initial baseline or the entirety of the [Hbtot] time course. Overall performance of each filter in the unmixed data is similar to that seen in the multispectral case (Fig. 2).

The αβ filter effectively captured transient spectral dynamics while rejecting noise in both tumor models (Fig. 5). This is in contrast to the SW’s lack of noise rejection, illustrated by the highly variable baseline signal, and the α’s lack of dynamic sensitivity, seen in the rounding-off of the gas challenge transition in the 700 nm channel (Fig. 5B, D: green line).

Fig. 5

Time courses of MSOT spine signal for representative channels in RENCA-luc mouse. A) Sum image of all channels at final time point, with ROIs used for analysis indicated. T: tumor, S: spine. B, C, D) Time courses of signal in selected channels based on respective filters. B) 700 nm, C) 800 nm, D) 850 nm. Red: sliding window filter, green: , blue: . The filter removed a substantial proportion of the noise from the dynamic signal while still accurately following the gas challenge. All channels are shown in Supp. Figures S1 (image comparisons) and S2 (time courses).

The dynamic performance of each filter for the RENCA-luc data is shown in Fig. 6, Fig. 7: Whereas the sliding window data are extremely noisy, preventing straightforward characterization of the response to gas challenge, the α filter clearly separated distinct subpopulation time courses within each ROI. However, the α filter introduced a substantial response time, taking minutes to reach steady state after each gas change, and consequently destroying quantitative information regarding the rate constant of the transition, as well as the transition time itself. In comparison, the αβ filter revealed conspicuous differences in both the rates and degree of response to gas challenge across both simulated and experimental data. These effects are more obvious in the transition region highlighted in Fig. 7: There was a dramatic decrease in the amount of noise in both the α and αβ filtered data, but only the filter preserved the rate and degree of response at each time point.

Fig. 6

Dynamic histogram images and selected histograms for each parameter-filter combination in RENCA-luc tumor-bearing mouse. ROIs defined in Fig. 5A. Upper panels: spine ROI. Lower panels: Tumor ROI. Each panel: Dynamic histogram images (left column) and corresponding histograms at selected time points (right column) for sliding window (top), (middle), and (bottom) filters. For all parameters, the use of either the or filters reduced the noise in the dynamic histograms, but the filter lagged behind across the entire ROI, as seen by the misalignment in the red and green histograms highlighted by the black arrows.

Fig. 7

Expanded time courses of air-oxygen transition in RENCA-luc tumor-bearing mouse. ROIs defined in Fig. 5A. Upper panels: Spine ROI. Lower panels: Tumor ROI. Left panels: [HbO2]MSOT. Right panels: [SO2]MSOT. A conspicuous reduction in noise is seen between the sliding window (each panel, top) and (each panel, bottom) time courses, allowing one to resolve the distinct evolution of different subpopulations within the ROI.

When applied to the XP373 tumor (Fig. 8, Fig. 9), the filter reveals heterogenous response to the gas challenge which is not visible in the sliding window data. Indeed, it is clear in the upper-left panel of Fig. 9 that there are both responsive and nonresponsive subpopulations of pixels within the spine ROI. The αβ filter also performed well on relatively constant data such as the [HbO2]MSOT channel in the tumor ROI, delineating a bimodal distribution which was obscured by noise in the SW filter. These general between-filter patterns are sustained when applying the filters to additional XP373 gas challenge datasets, as seen in Supplementary Figs. S3–5.

Fig. 8

Dynamic histogram images and selected histograms for each parameter-filter combination in XP373 tumor-bearing Mouse #1. ROIs defined in Supp. Fig. 3A. Upper panels: Spine ROI. Lower panels: Tumor ROI. For all derived parameters, the use of either the or filters reduced the amount of noise in the dynamic histograms, though the filter lagged behind across the entire ROI, as seen by the misalignment in the red and green histograms highlighted by the black arrows.

Fig. 9

Expanded time courses of air-oxygen transition in XP373 tumor-bearing mouse (Fig. 8). Upper panels: Spine ROI. Lower panels: Tumor ROI. Left panels: [HbO2]MSOT. Right panels: [SO2]MSOT. A conspicuous reduction in noise can be seen between the sliding window (each panel, top) and (each panel, bottom) time courses, allowing one to resolve the distinct evolution of different subpopulations within the ROI.

Both the and filters provided improved statistical performance in assessing whether an oxygen-induced change is significantly different from the air baseline, in both spatial and temporal dimensions (Fig. 10, Table 1). Areas where the signal to noise ratio is low are deemed as insignificant changes under analysis of the SW-filtered data, but both the and filters resolve those same areas as significant (Supplementary Fig. S6). In examining the transition over time from insignificant to significant, fitting of a logistic function to the binary classification time course yields an effective transition time. The analysis provided by the SW filter provided inconsistent fits for transition times, which were highly variable depending on the SNR of the local region. Though the filter provided a much sharper transition from non-significant to significant, it occurred at a substantial delay from the known switching time of 10.5 min. The filter yields consistent and sensible transition times between all ROIs.

Fig. 10

Significance testing of HbO2MSOT change due to gas challenge in RENCA-luc tumor-bearing mouse. Left: 4-pixel ROIs used. ROIs 1 and 2 are significant via the pixel-wise analysis for all filters, while ROI 3 is not significant for the SW filter (Table 1, Supplementary Fig. S2). Center: log10(p) values over time for each filter in each ROI. Right: Binary classification of change as significant (s.) or not (n.s.), as well as logistic curve fit and indicated switching time. The SW filter is conspicuously noisy throughout time for all ROIs, but particularly so in ROI 3, causing an unreasonable fit of the switching time. The filter denoises effectively and has a very sharp transition from non-significant to significant, but incurs a 0.5 min delay in the switching time. The filter denoises effectively while preserving the dynamic transition, giving a consistent switching time for all ROIs (see Table 1).

Table 1

Comparison of filter performance based on spatial and temporal analyses. The spatial analysis refers to the significance of the [HbO2]MSOT difference before and after the gas challenge transition within each ROI inFig. 10. The temporal analysis refers to the significance of the difference between a single [HbO2]MSOT image at a particular time point before the gas challenge transition and each successive [HbO2]MSOT image during the gas challenge transition.

ROI	Filter	Spatial analysis			Temporal analysis
		p-value	Scaled p-value	Significant	Mean p-value of t_1695:t_1700	tswitch (min)
1	SW	1.3×10^-8	7.8×10^-4	Yes	7.1×10^-6	11.53
	α	5.6×10^-24	3.5×10^-19	Yes	2.2×10^-6	11.95
	αβ	3.9×10^-13	2.4×10^-8	Yes	1.6×10^-6	11.41
2	SW	5.2×10^-9	3.2×10^-4	Yes	1.1×10^-6	11.77
	α	2.2×10^-26	1.4×10^-21	Yes	1.9×10^-6	12.07
	αβ	1.6×10^-14	9.9×10^-10	Yes	8.5×10^-7	11.66
3	SW	2.1×10^-6	0.13	No	5.2×10^-6	13.05
	α	1.6×10^-20	1.0×10^-15	Yes	9.3×10^-7	12.03
	αβ	9.4×10^-11	5.9×10^-6	Yes	3.5×10^-7	11.45

Discussion

The filter provides an effective balance between fidelity to spectral dynamics and noise rejection in a wide variety of contexts and under different spectral sampling patterns. Notably, the filter is robust to variations in the parameters and , and thus requires little tuning in order to achieve substantially improved results.

In scenarios where no appreciable change is expected, or where the dynamics are very slow relative to the sampling rate, the filter provides superior performance over the filter due to its effective noise-suppression. This may be used to accommodate piecewise-stationary statistics wherein large changes in the underlying signal occur in between periods of steady-state behavior, or when the noise power is large relative to the true signal variation.

Accommodating high-noise, low-variation environments may be accomplished by setting the measurement noise variance high relative to a lower process noise variance. This causes the filters to bias themselves towards previous estimates, providing benefits for such contexts as measuring hemodynamic cycling. Other dynamic imaging modalities such as real-time functional MRI [48,49] or calcium imaging [50] may benefit from the use of these filters as well, as experimental motifs such as neurofeedback require reduced data noise without incurring substantial analysis delays [51].

Due to the ability of the and filters to incorporate information from previous samples, it may be desirable to choose sampling patterns which acquire one image at each wavelength, rather than oversampling each wavelength sequentially. Although the latter case enables averaging to improve signal to noise ratio, it inherently degrades time resolution and can introduce substantial blurring due to motion [24]. Further, this work only explores the application of these filters to reconstructed images – it is probable that image quality itself could be improved by applying the methods to raw photoacoustic data, in order to better condition the reconstruction process itself, and we will investigate this matter in a forthcoming publication.

Motion may present a challenge for these filters due to their recursive structure. As the and filters are infinite impulse response (IIR) filters, an arbitrarily extreme observation or outlier will affect the output data at any later time. Rejection of spurious data could be accomplished by identifying frames which are misregistered and down-weighting the associated residuals. This identification process may be accomplished via a number of methods, ranging from simple differences of sequential frames for online analysis [24] to sophisticated clustering methods when post-processing is an option [[52], [53], [54]].

The work presented here assumes that , and are constant through time, as well as being homogeneous throughout the image and across channels. As measurement noise may vary between imaging channels and across a field of view, and different regions may have distinct dynamics, a single global tracking index may not be optimal. The filter performs well even with the use of global parameters, but it may nevertheless be desirable to enable these parameters to vary throughout the temporal, spatial, and spectral dimensions.

Although other methods such as Bayesian or Kalman filters might provide better quantitative and qualitative results, the filters described here are very efficient both in terms of computational complexity and required storage, and can be readily implemented in any processing or acquisition pipeline with minimal performance impact. Indeed, the SW and filters each require only a single copy of the spectral image, while the filter only requires the additional storage occupied by the image for substantially improved performance.

As these filters are extremely efficient, it is possible to simultaneously maintain several filters with different parameters, using the data likelihood for each in order to provide an improved estimate of the state [28,29]. This, in turn, may be used to modify the parameters for each filter, giving an efficient adaptive update scheme.

Table 2 shows the time complexities and storage requirements for each of the presented methods, as well as comparisons against the full Kalman filter (KF). We assume an overall multispectral image size of pixels over channels, a measured image size of N, and the use of global parameters for both and filters. For the KF, we assume non-sparse matrices and no control inputs. Timing for the various filters was performed using a single thread on an Intel Xeon E5-2680 v3 processor operating at 2.5 GHz. The algorithm-specific memory footprint was less than 5 MB for all filters tested.

Table 2

Comparison of filter performance for various filtering methods.

	Time Complexity	Storage Requirement	Time/frame (s)
SW	O(N)	N*M	8.41 (± 9.54) x 10−5
α	O(N)	N*M + 3	6.83 (±2.26) x 10−4
αβ	O(N)	2*N*M + 3	1.3 (± 0.42) x 10−−3
Kalman	O(N3M3)	4N2M2 + N2M + 2N2 +NM+2N	Not tested

The SW retains only the copy of the multispectral image, and overwrites the values in a given imaging channel with the measured values. It therefore has a time complexity of O(N) and a storage requirement of N*M. The filter retains a copy of the multispectral image as well as a single floating-point value to store the tracking index, or three values to represent the , , and values. The calculations are performed elementwise, and so the calculation has a time complexity of O(N), and a storage of N*M + 3. The filter retains a copy of the multispectral image as well as the dx image, giving it a storage requirement of 2*N*M + 3. The calculations still operate using basic arithmetic, and so the overall cost of an iteration is O(N).

The Kalman filter requires the storage of a N*M x N*M state-transition matrix F, a N*M x N*M covariance matrix P, a N*M x N*M process noise covariance matrix Q, measurement matrices H for each channel, yielding M * (N x N*M), a measurement noise covariance matrix R, of size N*N, the innovation covariance S, of size N*N, and must have storage for the Kalman gain matrix of size N*M x N. There must also be storage for the previous estimate, the measurement itself, and the innovation residuals, adding another M*N+2*N elements. The time complexity is dominated by matrix multiplications and the finding of a matrix inverse, and exceeds O(N3M3). Various optimizations may be performed in order to limit the complexity of the KF, but many induce their own overhead and do not lower the complexity to the linear time seen in these filters [55].

Conclusion

The and filters, originally developed for target tracking and now applied to spectral imaging, are able to effectively preserve the spatial, spectral, and temporal structure of MSOT data while improving the signal to noise ratio, and are able to do so in a highly efficient, online fashion. These filters are simple and easily implemented, theoretically optimal, and practically effective. The resultant increases in image quality enable new capabilities and better insight into highly dynamic biological systems, allowing improved hypothesis testing and decision making.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.