Feasibility of quantitative tissue characterization using novel parameters extracted from photoacoustic power spectrum

Biomedical Signal Processing and Control 57 (2020) 101719

Contents lists available at ScienceDirect

Biomedical Signal Processing and Control journal homepage: www.elsevier.com/locate/bspc

Feasibility of quantitative tissue characterization using novel parameters extracted from photoacoustic power spectrum Nikita Rathi a,∗ , Saugata Sinha a , Bhargava Chinni b , Vikram Dogra b , Navalgund Rao c a b c

Department of Electronics and Communication, Visvesvaraya National Institute of Technology, South Ambazari Road, Nagpur, 440010, India Department of Imaging Science, University of Rochester Medical Center, 601 Elmwood Ave, Box 648, Rochester, NY, 14642, USA Chester F. Carlson Centre for Imaging Science, Rochester Institute of Technology, 54 Lomb Memorial Drive, Rochester, NY, 14623, USA

a r t i c l e

i n f o

Article history: Received 2 April 2019 Received in revised form 19 September 2019 Accepted 10 October 2019 Keywords: Photoacoustic imaging Photoacoustic spectrum analysis Tissue characterization

a b s t r a c t Photoacoustic (PA) spectral parameters, obtained using conventional frequency domain analysis of photoacoustic power spectrum, bear complex relationship with absorber size, due to which it is difficult to recover quantitative information related to tissue microstructure accurately, using the conventional spectral parameters. In this study, we proposed to utilise two new spectral parameters namely number of lobes as well as average lobe width, which may change linearly with PA absorber radius, for recovering quantitative values of PA absorber size. Using an analytical model of PA power spectrum with bandlimited acquisition, we developed linear relationships between number of lobes and absorber radius as well as average lobe width and absorber radius for a wide range of radius values. These linear relationships were validated using simulations and experimentally acquired PA signal. The preliminary results of the simulation as well as phantom study showed that using the computed values of number of lobes as well as average lobe width in the linear relationships, absorber radius could be recovered with reasonable degree of accuracy. The proposed parameters were extracted from the multiwavelength PA data generated by freshly excised prostate tissue specimens from 30 human patients undergoing prostatectomy for biopsy confirmed prostate cancer. Statistical analysis showed that both the proposed parameters were significantly different between malignant and non malignant prostate as well as malignant and normal prostate. The preliminary findings showed that the proposed parameters can be used to differentiate diseased tissue pathology form normal tissue. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction Photoacoustic (PA) imaging, which is based on PA effect, is an emerging non invasive soft tissue imaging modality. In PA effect, light absorbing elements emit high frequency ultrasound (US) waves, known as PA waves, following absorption of nanosecond duration laser pulse. In PA imaging of soft tissue, the PA waves emitted by light absorbing tissue elements are captured by US transducers [1]. As it is a hybrid imaging modality between pure optical imaging and US imaging, it retains certain advantages of the both. Like pure optical imaging, it provides optical property based contrast which can be utilised to obtain functional information about different tissue pathologies. On the other hand, since its output is US waves, PA imaging does not suffer from rapidly degrading resolution with increasing depth inside soft tissue like

∗ Corresponding author. E-mail address: [email protected] (N. Rathi). https://doi.org/10.1016/j.bspc.2019.101719 1746-8094/© 2019 Elsevier Ltd. All rights reserved.

pure optical imaging. In other words, PA imaging overcomes the degrading resolution related limitation of pure optical imaging and mechanical contrast based poor tissue characterisation limitation of pure US imaging [1]. PA imaging systems acquire PA signals by scanning the tissue specimen or organ along different trajectories. Two dimensional grayscale PA images are reconstructed using the acquired PA signals. The PA pixels of the reconstructed grayscale PA images, represent PA signal strengths, generated at different tissue locations. The conventional approach of tissue characterisation uses PA image pixel values which depict the variation of optical absorption property inside tissue. However, the factors like light attenuation in overlying tissue, excessive scattering of light in soft tissue and frequency dependent attenuation of US waves inside soft tissue make PA pixel value a non linear function of optical absorption coefficient. Due to complex and non linear relationship between the PA signal and optical absorption property of tissue, researchers explore frequency domain analysis of PA signals for tissue characterisation in terms of size of PA absorbers. Similar to US imaging [2,3], in PA fre-

2

N. Rathi, S. Sinha, B. Chinni et al. / Biomedical Signal Processing and Control 57 (2020) 101719

quency domain analysis, three spectral parameters namely slope, midband fit and intercept are extracted from the calibrated power spectrum of the time dependent raw PA signal. Several research groups performed studies to utilise the spectral parameters for tissue characterisation. Using frequency domain analysis on ex vivo PA signal generated by ocular tissue, Silverman et al. [4] found that spectral slope and midband fit parameters were high around the iris pigment epithelium. Kumon et al. [5] differentiated between malignant and non-malignant prostate in an in vivo study of prostate cancer marine model by showing that midband fit and intercept parameters were significantly different between malignant and normal tissue regions. Yang et al. [6] showed that slope values corresponding to 2 different-diameter microsphere distributions were significantly different from each other using a phantom study. Relationship between PA spectral parameters with PA absorber radius and concentration was formulated by Xu et al. [7] using a simulation study. Wang et al. [8] showed analytical relationship between PA spectral parameters and size of the absorber. Using a phantom study, Chitnis et al. [9], demonstrated the correlation of slope, midband fit and effective absorber size extracted from the frequency spectrum with mean absorber diameter of the microspheres. Using ex vivo and in situ imaging studies on mouse models, Xu et al. [10] explored the feasibility of differentiating fatty liver from normal liver using spectral parameters extracted from PA power spectrum. Ex vivo studies with human patients suspected of suffering from prostate and thyroid cancer were performed by Sinha et al. [11,12] to differentiate between malignant, normal and benign tissue using PA spectral parameters. Relationship between vascular diameter of capillaries and spectral slope, extracted from PA power spectra were derived theoretically and later validated using PA imaging of phantoms by Gao et al. [13]. Using an analytical model, Xu et al. [14] derived the expression of PA power spectrum corresponding to an assembly of PA microspheres and investigated the relationship of absorber diameter with different spectral parameters for various concentrations of microspheres. Using simulation and experimental studies, Xu et al. [15] showed the PA spectral parameters were correlated with Gleason scores. Using simulations, Feng et al. [16] showed correlation between the MTT (mean trabecular thickness), an important factor determining the bone strength, with the spectral slope extracted from the PA power spectrum and validated their findings using an experimental study. Hysi et al. [17] showed that the correlation of PA spectral parameters and estimated oxygen saturation values can be used for monitoring treatment-induced changes in the tumour vasculature. Wang et al. [18] deduced the expression of PA power spectrum generated by a random mixture of particles with non-uniform sizes and showed that there is an approximate linear relation between the content of abnormal particles and the spectral slope. In a phantom study, they differentiated a trace of large particles mixed in small particles. Feng et al. [19] proposed the use of micro-ring ultrasonic detectors for improving the performance of photoacoustic spectral analysis to differentiate between different tissue pathologies by extending the signature of the spectral parameters from tissue level to cellular level. Gorey et al. [20] proposed to utilise the particular signature of PA spectrum to detect FVPTC tissue of a thyroid nodule. They also performed differentiation of other thyroid tissue variants using power spectra of PA signals obtained following empirical mode decomposition. So far, the primary focus of the research in the area of PA spectral parameters has been to establish and validate the fact that spectral parameters change significantly with the change of tissue microstructure. Only a few studies attempted to establish actual relationships, which could be used for computing values of parameters like absorber radius or absorber concentration, between different parameters of tissue microstructure and spectral parameters. As spectral slope is independent of optical absorption property of the PA absorber unlike the other two spectral parameters, it is

expected to be most useful parameter among the three parameters for detecting the change in tissue microstructure. Previous studies have shown that the relationship between spectral slope and absorber radius is non linear or piece wise linear i.e. linear for very small ranges of radius values. Determination of the exact non linear relationship using which the PA absorber radius values can be predicted accurately for given spectral slope values, is difficult. In this study, we introduced two new parameters namely number of lobes and average lobe width, which can be computed using power spectrum of PA signal, for accurate prediction of PA absorber radius values. The logarithmic power spectrum of PA signal generated by a spherical PA absorber possess periodic fluctuations which change as the radius of the PA absorber changes, as shown in Fig. 2. These fluctuations are referred here as lobes. Our hypothesis was both number of lobes as well as average lobe width, change linearly with the PA absorber radius. In this study, we attempted to verify our hypothesis by finding linear relationships between absorber radius and number of lobes as well as average lobe width. These linear relationships were validated using a simulation study as well as experimental data generated by spherical PA absorber. Finally, in a retrospective ex vivo PA imaging study with prostate tissue specimens from suspected prostate cancer patients, statistical analysis was performed to investigate whether the proposed parameters were significantly different between different prostate tissue categories. The rest of the paper is organized as follows. In Section 2, linear relationships between the proposed parameters and PA absorber radius were established considering bandlimited PA acquisition. Later in this Section, the linear relationships were verified using a simulation study. Section 3 contains the details of phantom and ex vivo prostate tissue studies. In Section 4, the results of simulation and experimental studies as well as the significance of the proposed parameters were discussed in detail, while Section 5 contains the final conclusion. 2. Analytical study In this section, the analytical expression of the PA power spectral density, generated by a single spherical PA absorber with uniform spherical absorption profile was shown considering a bandlimited US transducer. The expression of PA signal, generated by an assembly of spherical PA absorbers and captured by an US transducer, were given by Xu et al. [14]. Spatially, PA absorbers were distributed following a uniform discrete random distribution. Eq. (1) represents the expression of the time dependent PA signal, recorded by the US transducer. f (t) = [Nshape(t) ×





˛(ω, |r − r|) ]

∞   ∗

[U(r, , ˚, t)]D2 (, ˚)dd˚}g(r1 , r2 )dr

{ 0

(1)

− −

Here Nshape(t) represents the N shaped initial PA signal generated by each spherical PA absorber. Frequency dependent attenuation of the PA signal in soft tissue is represented by ␣, while U represents the discrete uniform random function. D2 and g represent the directivity associated with the receiving efficiency of the transducer and gating function. The N-shape profile is given by,

 v s

Nshape (t) = A −

Rs



t+1 , −

Rs

vs


Rs

vs

(2)

Where A is the amplitude of the N-shape profile, RS is the radius of PA absorber and vs is the speed of sound. Eq. (1) consists of two parts: (i) the source profile and (ii) the integration, which is

N. Rathi, S. Sinha, B. Chinni et al. / Biomedical Signal Processing and Control 57 (2020) 101719

3

As mentioned in [21], the spectrum of US transducer is mostly smooth, slow function with a maximum in the vicinity of the central frequency. We considered TF of US transducer as a Gaussian function. Combining Eq. (8) and (9), psd of acquired PA signal by a bandlimited US transducer can be formulated as,



Fig. 1. A LTI system.

sp determined only by the spatial factors including the spatial distribution of the PA sources, the spatial directivity function of the US transducer and the gating function. The integration part in Eq. (1) is defined as Spatial(r,t). Therefore, the power spectrum of recorded PA signal is,



∞

sp (ω) = [

Nshape(t) ∗ Spatial(r, t)e

−iωt

dt]

−∞



∞

conj[

Nshape(t) ∗ Spatial(r  , t)e−iωt dt]

(3)

−∞

sp (ω) =

 1 2

 1



2

Nshape (ω) conj[Nshape (ω)]



Spatial (r, ω) conj[Spatial (r, ω)]





∞

sp (ω) =

RNshape () e

−iω

(4)



∞

d

RSpatial () e

−∞

−iω

d

(5)

−∞

Where RNshape and RSpatial represent autocorrelations of Nshape(t) and Spatial (r,t). From Eq. (2) we have +

2Rs

s

A2 (−

=

RNshape ()

s s t + 1)[− (t − ) + 1]dt Rs Rs

−

2Rs <<0 s

(6)

0

2Rs s A2 (−

=

s s t + 1)[− (t − ) + 1]dt Rs Rs

0≤<

2Rs s



The contribution of the source profile to the power spectrum is thereby,



sp

Nshape

∞

(ω) =

RNshape () e−iω d

(7)

−∞

 sp

  2

2

Nshape

cos

(ω) = A

 2ωR  s

vs

2

ω4 Rvss −

4 ω3 Rvss

sin

2 + 2 ω

 2ωR  s

vs

+

 2 2 − 2 ω2 ω4 Rvss (8)

Eq. (8) corresponds to the power spectrum of initial PA signal generated by a PA absorber. To ensure penetration at significant depth inside soft tissue, US as well as PA acquisition systems are generally bandlimited because of frequency dependent attenuation. To incorporate the effect of bandlimited acquisition, we assumed the PA acquisition system as a linear time invariant (LTI) system with bandlimited transfer function, as shown in Fig. 1. Eq. (9) provides the output power spectral density (psd) of the bandlimited PA acquisition system, where SX (w), SY (w) and H(w) are the input psd, output psd and TF of the system.



2

SY (ω) = H (ω) SX (ω)

(9)

out (ω)

cos

= A2

 2ωR  s

vs

 

2

2

ω4 Rvss −

4 ω3 Rvss

sin

2 + 2 ω

+



2 2 − Rs 2 ω2 4 ω vs

 2ωR   s

vs

e

−(ω−m)2 2 2

2

(10)

Where m and ␴2 represent mean and variance of gaussian distribution. The values of m and ␴2 depend upon the central frequency and bandwidth of the transducer. As we wanted to validate the findings of the analytical power spectrum model using simulations as well as experimental PA signals, the TF of the US transducer used in the analytical model was required to be similar to the US transducer used in our experimental PA imaging system. The US transducer, used in our PA imaging setup, had 5 MHz usable bandwidth [10 dB] (2.4 MHz–7.4 MHz) and 5 MHz centre frequency. We selected m and ␴ as 5 and 1.65 respectively and power spectrum corresponding to a single PA absorber was obtained using Eq. (10) considering bandlimited acquisition. As mentioned before, logarithm of PA power spectrum generates periodic fluctuations which we refer here as lobes. Fig. 2a shows the log PA power spectrum corresponding to spherical PA absorber of 0.5 mm radius (blue line) and 1 mm radius (red line), computed using Eq. (10). As it can be seen in Fig. 2a, with the change of absorber radius, both number of lobes as well as width of each lobe changed. In this study, our objective was to compute the PA absorber radius using the values of number of lobes and average lobe width, extracted from log power spectrum of PA signal generated by a spherical absorber. As seen in Fig. 2a, each lobe in the PA log power spectrum is defined by extreme points i.e. two minima and one maximum point. Based on this observation, we decided to compute number of lobes in terms of number of extreme points and lobe width in terms of distance between alternate extreme points. For locating the extreme points, we decided to find out the stationary points in the differentiated log power spectrum. Once the logarithm of power spectrum of PA signal was computed, it was differentiated to locate the extreme points. For differentiation, we computed detail wavelet coefficients using db1 mother wavelet. Fig. 2b shows the detail spectra of the logarithms of two different PA power spectra, shown in Fig. 2a. The zero crossing points in the detail spectrum were identified as the extreme points in the log PA power spectrum. Number of lobes was computed using the number of extreme points in the usable bandwidth while average lobe width was calculated using the distances between alternate extreme points in the usable bandwidth. Both number of lobes as well as average lobe width were computed for eight radius values starting from 0.25 mm to 2 mm at an interval of 0.25 mm. Fig. 3a and 3b show the plots of number of lobes against different radius values and average lobe width against different radius values respectively. The figures show that while number of lobes increased linearly with PA absorber radius, average lobe width decreased linearly with PA absorber radius. Eq. (11) and (12) are the Equations of the linear regression models, fitted to the points in Fig. 3a and b. y = 14.38Rs + 2.071

(11)

z = −0.415Rs + 1.123

(12)

Where, y and z are number of lobes and average lobe width respectively.

4

N. Rathi, S. Sinha, B. Chinni et al. / Biomedical Signal Processing and Control 57 (2020) 101719

Fig. 2. (a) Log power spectra corresponding to spherical PA absorbers of 0.5 mm radius (blue line) and 1 mm radius (red line), (b) Detail wavelet spectra corresponding to (a) in the usable bandwidth.

Fig. 3. (a) Number of lobes vs radius values of PA absorber, (b) Average lobe width vs radius values of PA absorber.

Fig. 4. Spectral slope vs radius values of PA absorber.

To compare the relation of spectral slope with the PA absorber radius with that between number of lobes and absorber radius as well as average lobe width and absorber radius, we also computed spectral slope values for the same range of radius values. We computed PA power spectrum, calibrated for the bandlimited acquisition, using Eq. (8). Spectral slope was determined by finding the slope of the linear regression model best fitted to the portion of log power spectrum with usable bandwidth. In Fig. 4, different values of spectral slope were plotted against the corresponding values of absorber radius. It can be seen that the slope values changed very slowly with radius values for the given range. Eq. (13) shows the linear regression model fitted to the points in Fig. 4. w = −.007Rs − .0174

(13)

Where, w is spectral slope. To verify the linear relationships between the proposed parameters and the PA absorber radius, we performed a small simulation study. Time domain PA signals corresponding to spherical PA absorbers with particular radii were simulated using the expression provided by Diebold et al. [22]. The same eight radius values covering the same range as used to generate Fig. 3a and Fig. 3b,

were considered here. The simulated PA signals were convolved with the impulse response of a 32 element 1D US transducer array (centre frequency, 5 MHz; bandwidth, 60%; pitch, 0.7 mm; elevation, 1 mm; Olympus NDT, State College, PA), which was used in our experimental PA imaging setup, to simulate the bandlimited PA signals. Log power spectrum of each bandlimited PA signal was computed and the two proposed parameters were calculated using the same process, used for computing the parameters from the analytical model of bandlimited PA logarithmic power spectrum. Fig. 5 shows the flowchart for computation of absorber radius from the proposed parameters. Fig. 6 shows the plots of computed values of the parameters for the simulated PA Aline signals vs absorber radius and the predicted values of the parameters, obtained using Eq. (11) and (12) vs absorber radius. The average error between the parameter values obtained for simulated PA signals and the parameter values obtained using Eq. (11) and (12) are −0.252 and −.045. 3. Experimental study 3.1. Experimental validation with phantom data To correlate the findings of the analytical as well as simulation studies, we computed the same parameters using the PA signal generated by a spherical phantom with known radius. Using the linear relationships between the parameters and the PA absorber radius, developed in the analytical study, we calculated the radius of the PA absorber and compared the computed radius value with the actual radius value. We used a mustard seed as a spherical PA absorber with radius 0.65 mm. Using the acoustic lens-based PA imaging system, developed in our laboratory, we acquired the PA signal generated by the mustard seed. Although the details of the PA imaging setup have been mentioned in our earlier reports [11,12], for the purpose of continuity, a short and concise description of the PA imaging system is provided here. The PA imaging system, as shown in Fig. 7

N. Rathi, S. Sinha, B. Chinni et al. / Biomedical Signal Processing and Control 57 (2020) 101719

5

analytical expression (blue line) and experimentation (red line). The number of lobes and average lobe width within the usable bandwidth range were 10 and 1.08 MHz respectively. These values were plugged into Eq. (11) and Eq. (12) to calculate the PA absorber radius. The PA absorber radius values, calculated using the number of lobe and average lobe width were 0.55 mm and 0.58 mm, while actually the PA absorber radius was 0.65 mm. The value of the spectral slope, which was computed using the log power spectrum of the PA A line signal, calibrated with the log power spectrum of transducer array in the usable bandwidth range, was −0.019 dB/MHz. Putting the value of the spectral slope in Eq. (13), the radius of the mustard seed was found out to be 0.23 mm.

3.2. Experiment with human prostate data

Fig. 5. Flowchart for computation of radius of PA absorber.

was developed following a transillumination geometry in which the laser light source and PA detection system were placed on the opposite sides of the sample holder. The laser light source was the fibre optic bundle, connected with the multiwavelength nanosecond pulse laser, while the PA detection system was a water filled cylinder containing a single element acoustic lens and an US transducer array. The PA waves generated by the sample following laser light irradiation, were focussed by the acoustic lens on the transducer array. The transducer array was connected to a custom designed data acquisition system which amplified, digitised and stored the PA signal in the computer. The sample holder containing the mustard seed, embedded in a Polyvinyl chloride plastisol pad, was filled with water before PA imaging. Once we acquired the radio frequency PA signals generated by the mustard seed, we selected the A line PA signal corresponding the centre element of the transducer array as the centre element was placed directly on top of the mustard seed. Fig. 8 shows 2D image of a mustard seed, 2D PSF of mustard seed, 1D Aline PA signal acquired from mustard seed and corresponding log magnitude spectrum. We computed log power spectrum of the A line PA signal and calculated number of lobes as well as the average lobe width in the usable bandwidth range. For identifying the lobes in the log power spectrum, we computed its detail wavelet coefficients using db1 as mother wavelet, in the same way as mentioned in the analytical study. Fig. 9 shows log power spectrum of acquired PA signal with 0.65 mm radius obtained by simulation (black line),

We extracted the parameters from multiwavelength PA data acquired in an ex vivo PA imaging study with actual human patients undergoing prostatectomy for biopsy confirmed prostate cancer. The PA imaging protocol, followed in our laboratory for ex vivo imaging of human patients, was mentioned in detail in our earlier reports [11,12]. In this report, a brief summary of the imaging protocol is mentioned for the sake of continuity. Following surgery, the surgical pathology laboratory received the exercised prostate gland, which was marked with ink and divided into 2–5 mm thick sections, by a genitourinary pathologist. PA data at multiple wavelengths, generated by an excised tissue specimen with glossy visible nodule, were acquired in our laboratory. After PA imaging was completed histopathology was performed on the excised prostate specimens. The complete procedure up to histopathology of the prostate specimen was finished within 1 hour after the surgery. The complete ex vivo PA imaging protocol was approved by the genitourinary pathologist. During PA imaging, the laser intensity on the tissue specimen was maintained around 5 mJ/cm2 , well below ANSI guideline for safe human exposer. The different wavelengths were chosen in such a manner that each wavelength coincided with the peak of the absorption spectrum of a particular chromophore, expected to be present in the tissue. Originally, the ex vivo study was performed at 5 different wavelengths, while in our retrospective study, PA data acquired at 3 different wavelengths were considered. Among the 3 different wavelengths, 760 nm and 850 nm are the absorption spectra peaks of deoxyhaemoglobin and oxyhaemoglobin with near infra-red (NIR) region, while at 800 nm, specific absorption coefficient of both deoxyhaemoglobin and oxyhaemoglobin are equal. Our study involved thirty patients with 53 ROIs. According to the histopathology analysis performed by the genitourinary pathologist, among the 53 ROIs, 19 were malignant, 8 were BPH, and 26 were normal. Here, we defined one more tissue category as non-malignant and under this category all ROIs corresponding to normal and BPH tissue were combined. Using the parameter values, we performed a 2-sample, 2-tailed t test to find out whether these parameters were statistically significantly different corresponding to ROIs of different tissue types. Fig. 10 shows the results of the t test performed at a 5% significance level with the values of the 2 parameters corresponding to malignant versus normal, malignant versus BPH, BPH versus normal, and malignant versus non-malignant. From Fig. 10, it can be seen that the 2 parameters were significantly different (P < .05) corresponding to the malignant and normal prostate tissue at the 760-nm, 800-nm and 850-nm wavelengths. Between malignant prostate and BPH tissue, number of lobes and average lobe width parameters were significantly different (P < .05) at all 3 wavelengths. For BPH versus normal prostate, none of the parameters, except number of lobes at 850 nm, were significantly different (P < .05). Both the parameters at all 3

6

N. Rathi, S. Sinha, B. Chinni et al. / Biomedical Signal Processing and Control 57 (2020) 101719

Fig. 6. (a) Number of lobes vs radius values of PA absorber, (b) Average lobe width vs radius values of PA absorber.

Fig. 7. Experimental set up.

wavelengths were significantly different (P < .05) corresponding to malignant and non-malignant prostate.

4. Discussion In this study, two new parameters were introduced for quantitative measurement of size of spherical PA absorbers. Using the analytical model of PA power spectrum provided by Xu et al. [14], we computed the number of lobes as well as average lobe width corresponding to spherical PA absorbers with different radius values ranging from 0.25 mm to 2 mm. In the analytical model of PA power spectrum, we introduced bandlimited PA acquisition. For computing the values of number of lobes as well as average lobe width, we considered the acquisition bandwidth to be 2.4 MHz to 7.4 MHz, which was same as that of the US transducer array in our actual PA imaging system. Fig. 3a and b show that number of lobes increased while average lobe width decreased with increasing PA absorber radius values. As it was observed that both number of lobes and average lobe width changed almost linearly with the change in PA absorber radius value, linear regressions were performed to find out relationships between PA absorber radius and number of lobes as well as average lobe width, as shown in Eq.

(11) and (12). Using the analytical model of PA power spectrum, we computed the values of spectral slope for the same set of spherical PA absorber radius values, used earlier for computing number of lobe and average lobe width computation. Fig. 4 shows that the spectral slope values remained almost unchanged for different PA absorber radius values which was further corroborated by the very small slope of the linear equation between PA absorber radius and spectral slope value, obtained using linear regression. The relationships obtained using the analytical PA power spectrum model were verified using a simulation study where time domain PA signals, generated by spherical PA absorbers of different radius values and acquired by a bandlimited US transducer, same as that was used in our experimental PA imaging system, were simulated. The number of lobes and average lobe width were calculated from the logarithmic power spectra of the simulated PA signals. Fig. 6a and b show the plots of number of lobes vs PA absorber radius as well as average lobe width vs PA absorber radius, calculated using the simulated PA signals as well as the analytical model. The average error between the parameter values, computed using the analytical PA power spectrum model and simulated PA signals are −0.25 for number of lobes and −0.045 for average lobe width.

N. Rathi, S. Sinha, B. Chinni et al. / Biomedical Signal Processing and Control 57 (2020) 101719

7

Fig. 8. (a) 2D image of a mustard seed, (b) 2D PSF of mustard seed, (c) 1D PA signal acquired from mustard seed, (d) log magnitude spectrum of (c).

For experimental verification of the relationships between PA absorber radius and number of lobes as well as average lobe width, we acquired PA signal, generated by a mustard seed using our PA imaging system. The values of the two proposed parameters as well as the spectral slope were computed from the logarithmic power spectrum of the PA signal. Using the computed values of these three parameters in Eqs. (11), (12) and (13), radius of the mustard seed was recovered. The radius values of the mustard seed, extracted using the computed values of number of lobes, average lobe width and spectral slope were 0.55 mm, 0.58 mm and 0.23 mm, while the actual radius of the mustard seed was 0.65 mm. The percentage errors, related to the recovery of the mustard seed radius using number of lobes, average lobe width and spectral slope were 15.38%, 10.77% and 64.62%. These results show that while it was possible to recover absorber radius using number of lobes and average lobe width with reasonable accuracy, performance of spectral slope was quite poor for predicting the PA absorber radius. Multiwavelength PA data, generated by freshly excised prostate tissue specimens from human patients, were analysed to test whether the proposed parameters were significantly different between different tissue pathologies. Two sampled two tailed student t tests were performed to check whether the number of lobes and average lobe width were significantly different between malignant, benign, normal and non-malignant prostate tissue. As mentioned earlier, PA data at three different wavelengths were analysed here. So, statistical analyses were done for a total of six different parameters i.e. number of lobes and average lobe width at three different wavelengths. These six parameters were tested for four pairs of different prostate tissue types, namely malignant vs normal, malignant vs BPH, BPH vs normal, and malignant vs nonmalignant. For malignant versus normal prostate tissue, malignant versus BPH and malignant vs non-malignant, 6 out of 6 parameters were significantly different (P < .05), whereas for BPH versus normal prostate, none of the parameters were significantly different (P < .05) except number of lobes at 850 nm. As far as the number

Fig. 9. Log power spectrum of acquired PA signal with 0.65 mm radius obtained by simulation (black line), analytical expression (blue line) and experimentation (red line).

of significantly different parameters are concerned, it is evident that the proposed parameters were least effective in differentiating BPH from normal prostate. Although, all six parameters were significantly different for the other three tissue category pairs, the p values corresponding to malignant vs non-malignant prostate were least. For differentiating different tissue pathologies, performance of number of lobes was better than that of average lobe width, as barring once, number of lobes provided lower p values than that of average lobe width consistently for all pairs of tissue categories at three wavelengths. Comparing among different wavelengths, it can be seen from Fig. 10 that both average lobe width and number of lobes performed best at 850 nm. In the conventional frequency analysis of PA signals, three parameters namely slope, midband fit and intercept of the straight line, best fitted to the portion of calibrated power spectrum in the usable bandwidth region, are computed. Among these three parameters, both midband fit as well as intercept depend on the size as well as the optical absorption of the PA absorber. As spectral slope does not depend on optical absorption property, it is generally con-

8

N. Rathi, S. Sinha, B. Chinni et al. / Biomedical Signal Processing and Control 57 (2020) 101719

Fig. 10. Results of the t test performed at a 5% significance level with the values of the 2 parameters corresponding to (a)malignant versus normal, (b) malignant versus BPH, (c) BPH versus normal and (d) malignant versus non-malignant.

sidered as the primary feature to detect any change in PA absorber size. Both of the proposed parameters do not depend on the optical absorption property i.e. even if the optical property changes as far as PA absorber size and acquisition system remain unchanged, the parameter values will be unchanged. Findings of our simulation as well as phantom study show that both of the proposed parameters performed much better than the spectral slope for predicting the PA absorber radius. This might have happened due to the bandlimited nature of the US transducer we considered. In preliminary simulation, we observed that spectral slope changes with PA absorber radius when the acquisition bandwidth extends to significantly low frequency limit. In other words, the effect of PA absorber radius change will be reflected in the spectral slope value when the usable bandwidth contains a significant portion of the primary or first lobe of the PA power spectra. In our case, the lower limit of the usable bandwidth was 2.4 MHz, which did not include significant portion of the first lobe PA power spectrum, as seen in the Fig. 2. On the other hand, the proposed parameters will change with the change of absorber radius as long as the acquisition bandwidth contains multiple lobes. Results of the earlier experimental and simulation studies, performed by different research groups showed that the relationship between spectral slope and radius of the regular shape PA absorber was nonlinear or piece wise linear. Formulation of the exact nonlinear relationship using which accurate quantitative recovery of the PA absorber radius value from measured spectral slope is quite difficult. On the other hand, as shown in Fig. 3a and Fig. 3b, the proposed parameters have simple linear relationships with absorber radius, using which absorber radius can be recovered from computed values of the proposed parameters with reasonable degree of accuracies. Spectral slope values are extracted from PA power spectrum which needs to be calibrated using the power spectrum of the system transfer function. In case of the proposed parameters, calibration does not change the parameter values as the number of lobes as well as average lobe width remain unchanged after calibration. Fig.11 shows the calibrated and uncalibrated logarithmic

Fig. 11. Calibrated (red line) and uncalibrated (blue line) logarithmic power spectra of the PA signal generated by a spherical PA absorber with radius 1 mm.

power spectra of the PA signal generated by a spherical PA absorber with radius 1 mm. As the number of lobes and width of each lobe do not change after calibration, values of both number of lobes and average lobe width are not changed following calibration. Unlike spectral slope computation, calibration is not required for computing the values of the proposed parameters. 5. Conclusion In this study, we investigated the feasibility of performing quantitative tissue characterisation using number of lobes and average lobe width, computed from PA power spectrum. The preliminary results showed that using the linear relationships between these two parameters and PA absorber radius, established here, the PA absorber radius could be computed with reasonable accuracy from the measured values of these parameters. The parameter values, extracted from the power spectrum of multiwavelength PA data generated by freshly excised prostate specimens from human patients, were significantly different between malignant and non malignant prostate tissue as well as malignant and normal prostate tissue.

N. Rathi, S. Sinha, B. Chinni et al. / Biomedical Signal Processing and Control 57 (2020) 101719

Acknowledgment This work is an outcome of the R&D work undertaken in the project under the Visvesvaraya PhD Scheme of Ministry of Electronics & Information Technology, Goveronment of India, being implemented by Digital India Corporation (formerly Media Lab Asia). Declaration of Competing Interest We confirm that this work is original and has not been published elsewhere nor is it currently under consideration for publication elsewhere. References [1] P. Beard, Biomedical photoacoustic imaging, Interface Focus 1 (4) (2011) 602–631, http://dx.doi.org/10.1098/rsfs.2011.0028. [2] F.L. Lizzi, E.J. Feleppa, S.K. Alam, C.X. Deng, Ultrasonic spectrum analysis for tissue evaluation, Pattern Recognit. Lett. 24 (4-5) (2003) 637–658, http://dx. doi.org/10.1016/S0167-8655(02)00172-1. [3] F.L. Lizzi, M. Greenebaum, E.J. Feleppa, M. Elbaum, D.J. Coleman, Theoretical framework for spectrum analysis in ultrasonic tissue characterization, J. Acoust. Soc. Am. 73 (4) (1983) 1366–1373, http://dx.doi.org/10.1121/1. 389241. [4] R.H. Silverman, F. Kong, H.O. Lloyd, Y.C. Chen, Fine-resolution photoacoustic imaging of the eye Photons Plus Ultrasound: Imaging and Sensing 2010, vol. 7564, International Society for Optics and Photonics, 2010, http://dx.doi.org/ 10.1117/12.842471, February. p. 75640Y. [5] R.E. Kumon, C.X. Deng, X. Wang, Frequency-domain analysis of photoacoustic imaging data from prostate adenocarcinoma tumors in a murine model, Ultrasound Med. Biol. 37 (5) (2011) 834–839, http://dx.doi.org/10.1016/j. ultrasmedbio.2011.01.012. [6] Y. Yang, S. Wang, C. Tao, X. Wang, X. Liu, Photoacoustic tomography of tissue subwavelength microstructure with a narrowband and low frequency system, Appl. Phys. Lett. 101 (3) (2012), 034105, http://dx.doi.org/10.1063/1.4736994. [7] G. Xu, I.A. Dar, C. Tao, X. Liu, C.X. Deng, X. Wang, Photoacoustic spectrum analysis for microstructure characterization in biological tissue: a feasibility study, Appl. Phys. Lett. 101 (22) (2012), 221102, http://dx.doi.org/10.1063/1. 4768703. [8] S. Wang, C. Tao, Y. Yang, X. Wang, X. Liu, Theoretical and experimental study of spectral characteristics of the photoacoustic signal from stochastically distributed particles, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 62 (7) (2015) 1245–1255, http://dx.doi.org/10.1109/TUFFC.2014.006806. [9] P.V. Chitnis, J. Mamou, A. Sampathkumar, E.J. Feleppa, Spectrum analysis of photoacoustic signals for tissue classification Photons Plus Ultrasound: Imaging and Sensing 2014, vol. 8943, International Society for Optics and Photonics, 2014, http://dx.doi.org/10.1117/12.2037929, March. p. 89432J.

9

[10] G. Xu, Z. Meng, J. Lin, P. Carson, X. Wang, Functional pitch of a liver: fatty liver disease diagnosis with photoacoustic spectrum analysis Photons Plus Ultrasound: Imaging and Sensing 2014, vol. 8943, International Society for Optics and Photonics, 2014, http://dx.doi.org/10.1117/12.2037384, March. p. 89431G. [11] S. Sinha, N.A. Rao, B.K. Chinni, V.S. Dogra, Evaluation of frequency domain analysis of a multiwavelength photoacoustic signal for differentiating malignant from benign and normal prostates: ex vivo study with human prostates, J. Ultrasound Med. 35 (10) (2016) 2165–2177, http://dx.doi.org/10. 7863/ultra.15.09059. [12] S. Sinha, N. Rao, B. Chinni, J. Moalem, E.J. Giampolli, V. Dogra, Differentiation between malignant and normal human thyroid tissue using frequency analysis of multispectral Photoacoustic images, in: 2013 IEEE Western New York Image Processing Workshop (WNYIPW), IEEE, 2013, pp. 5–8, http://dx. doi.org/10.1109/WNYIPW.2013.6890979, November. [13] X. Gao, C. Tao, X. Wang, X. Liu, Quantitative imaging of microvasculature in deep tissue with a spectrum-based photo-acoustic microscopy, Optics Lett. 40 (6) (2015) 970–973, http://dx.doi.org/10.1364/OL.40.000970. [14] G. Xu, J.B. Fowlkes, C. Tao, X. Liu, X. Wang, Photoacoustic spectrum analysis for microstructure characterization in biological tissue: analytical model, Ultrasound Med. Biol. 41 (5) (2015) 1473–1480, http://dx.doi.org/10.1016/j. ultrasmedbio.2015.01.010. [15] G. Xu, M.C. Davis, J. Siddiqui, S.A. Tomlins, S. Huang, L.P. Kunju, J.T. Wei, X. Wang, Quantifying Gleason scores with photoacoustic spectral analysis: feasibility study with human tissues, Biomed. Optics Express 6 (12) (2015) 4781–4789, http://dx.doi.org/10.1364/BOE.6.004781. [16] T. Feng, J.E. Perosky, K.M. Kozloff, G. Xu, Q. Cheng, S. Du, J. Yuan, C.X. Deng, X. Wang, Characterization of bone microstructure using photoacoustic spectrum analysis, Optics Express 23 (19) (2015) 25217–25224, http://dx.doi.org/10. 1364/OE.23.025217. [17] E. Hysi, L.A. Wirtzfeld, J.P. May, E. Undzys, S.D. Li, M.C. Kolios, Photoacoustic signal characterization of cancer treatment response: Correlation with changes in tumor oxygenation, Photoacoustics 5 (2017) 25–35, http://dx.doi. org/10.1016/j.pacs.2017.03.003. [18] S. Wang, C. Tao, X. Gao, X. Wang, X. Liu, Quantitative photoacoustic examination of abnormal particles hidden in a mixture of particles with non-uniform sizes, Optics Express 23 (25) (2015) 32253–32260, http://dx.doi. org/10.1364/OE.23.032253. [19] T. Feng, Q. Li, C. Zhang, G. Xu, L.J. Guo, J. Yuan, X. Wang, Characterizing cellular morphology by photoacoustic spectrum analysis with an ultra-broadband optical ultrasonic detector, Optics Express 24 (17) (2016) 19853–19862, http://dx.doi.org/10.1364/OE.24.019853. [20] A. Gorey, P.M. Jacob, D.T. Abraham, R. John, M.T. Manipadam, M.S. Ansari, G.C. Chen, S. Vasudevan, Differentiation of malignant from benign thyroid nodules using photoacoustic spectral response: a preclinical study, Biomed. Phys. Eng. Express 5 (3) (2019), 035017 https://orcid.org/0000-0002-3092-4621. [21] O. Michailovich, D. Adam, Shift-invariant, DWT-based “projection” method for estimation of ultrasound pulse power spectrum, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 49 (8) (2002) 1060–1072, http://dx.doi.org/10.1109/ TUFFC.2002.1026018. [22] G.J. Diebold, T. Sun, Properties of photoacoustic waves in one, two, and three dimensions, Acta Acust. United Acust. 80 (4) (1994) 339–351.