## 1.

## Introduction

Speckles, a high-contrast, fine-scale, granular image texture, prevail in all scattering-based coherent imaging modalities, including laser imagery,^{1} ultrasonography (US),^{2}^{,}^{3} synthetic aperture radar,^{2} and optical coherence tomography.^{3} When multiple scatterers are randomly distributed in a resolution cell of the imaging system, they cannot be resolved. Consequently, the partial waves scattered from individual scatterers interfere and form a speckle. Speckle grains create image artifacts that do not correspond to the real structure of the object.^{1} Many efforts have been made to mitigate speckle artifacts at the expense of system complexity, imaging time, or spatial resolution.

Speckles become fully developed when coherent partial waves have phases completely randomized within a full range of $2\pi $. Here, coherence refers to correlation instead of in-phase. For wideband waves, the phase is associated with the center frequency. Speckles, however, are under-developed in some circumstances. For example, if mirror surfaces are imaged by optical coherence tomography, or smooth bone surfaces are imaged by ultrasonography, we observe a phenomenon referred to as ‘specular reflection.’^{4} Specular reflection is formed by constructive interference of coherent partial waves that have phases randomized only within $[0,\pi ]$ or less. For scattering-based coherent imaging, two conditions must be met to ensure specular reflection. First, all scatterers on the surface must have similar properties so that the polarity changes to the incident wave due to backscattering are the same. Second, the boundary roughness must be less than $\lambda /4$—“with $\lambda $ being the center wavelength”—so that the phase-delay variations due to scatterer spatial distribution differences are within $\pi $.

Photoacoustic (PA) tomography (PAT) is a promising imaging technology that exploits optical absorption contrast from absorbers, such as hemoglobin,^{5}^{–}^{10} melanin,^{6} water,^{11} organic dyes,^{12} and nanoparticles,^{13} using ultrasonic detection.^{14} In PAT, tissues are usually irradiated by a short-pulsed laser. Absorbed light is converted into heat, which is further converted to a positive pressure rise via thermoelastic expansion. The initial pressure rise then propagates as ultrasonic waves, which are detected by one or more ultrasonic transducers. Then, all of the received ultrasonic signals are used to form an image.

As a coherent imaging modality based on optical absorption, PAT is generally devoid of speckles. For the following reasons identified in our previous work,^{15} speckles in PAT are usually suppressed by prominent boundary buildups, referred to as “specular emission” here, via a mechanism similar to that of specular reflection. First, because the initial PA pressure rise is always positive, the propagated PA wave from any finite-sized optical absorber starts with a positive pressure and ends with a negative pressure.^{16} In addition, because most PAT systems are linear and shift-invariant, the changes to the polarities of these partial waves due to imaging systems are the same. As a result, these partial waves possess the same initial polarities. Therefore, PAT naturally satisfies the first condition for “specular reflection” except that “specular emission” is a more accurate description.

Second, the optical absorbing targets in biological tissues, such as blood vessels, usually have smooth surfaces on the scale of the wavelengths of megahertz ultrasonic waves. Thus, the PA partial waves from most surface absorbers are generated within a region of $\lambda /2$ in thickness at the boundary of the absorbing target, yielding phase-delay variations within $\pi $. Consequently, the second condition for specular emission is also satisfied, and the constructive interference among the partial waves from individual absorbers leads to boundary buildups. By contrast, the PA partial waves generated from the interior of a sufficiently large absorbing target have phase-delays randomized over a full range of $2\pi $. Therefore, the second condition for specular emission is violated, and speckles are formed from the interior of the absorbing targets. We also found that the interior speckles are usually suppressed by the prominent boundary buildups. In a special case, if the absorbing object is smaller than $\lambda /2$, then the phase-delay variations of the PA partial waves from individual absorbers are smaller than $\pi $. Therefore, the object appears as a solid area in PA images, as a result of the constructive interferences. However, a natural question is whether the boundary signals in PAT can become fully developed speckles if the absorbing target has sufficiently rough boundaries, which is addressed in this paper.

Previously,^{15} we used the same linear model to study speckles for both PAT and US. An optically absorbing or acoustically scattering target is modeled as a collection of statistically independently and randomly distributed unit optical absorbers or scatterers, which are much smaller than the imaging wavelength. The particles in a volume $V$ are distributed with a number density of $\rho $. The expected instantaneous power of PA or US signals from the particles within the volume $V$ is given^{15} by:

## Eq. (1)

$$\langle P(t)\rangle =({\overline{C}}^{2}+{\sigma}_{C}^{2})\rho {\int}_{V}{[\tilde{h}(\stackrel{\rightharpoonup}{r},t-|\stackrel{\rightharpoonup}{r}|/c)]}^{2}d{\stackrel{\rightharpoonup}{r}}^{3}+{\overline{C}}^{2}{\rho}^{2}{[{\int}_{V}\tilde{h}(\stackrel{\rightharpoonup}{r},t-|\stackrel{\rightharpoonup}{r}|/c)d{\stackrel{\rightharpoonup}{r}}^{3}]}^{2}\mathrm{.}$$The first term on the right-side of Eq. (1) is the sum of the powers of the partial waves generated from all of the particles. This uncorrelated contribution to the total power represents the random fluctuations—speckles. The second term is the correlated contribution to the total power, which is responsible for the specular emission in PAT or the specular reflection in US. The specular emission can occur only at the boundaries, because ${\int}_{V}\tilde{h}(\stackrel{\rightharpoonup}{r},t-|\stackrel{\rightharpoonup}{r}|/c)d{\stackrel{\rightharpoonup}{r}}^{3}$ is always zero inside the absorbing target larger than the resolution cell.^{15} The speckle visibility in PAT is defined as the ratio of the square root of the average speckle power to the average magnitude of boundary features:

## Eq. (2)

$$V=\sqrt{\frac{({\overline{C}}^{2}+{\sigma}_{C}^{2})\rho {\int}_{V}{[\tilde{h}(\stackrel{\rightharpoonup}{r},t-|\stackrel{\rightharpoonup}{r}|/c)]}^{2}d{\stackrel{\rightharpoonup}{r}}^{3}}{{\overline{C}}^{2}{\rho}^{2}{[{\int}_{V}\tilde{h}(\stackrel{\rightharpoonup}{r},t-|\stackrel{\rightharpoonup}{r}|/c)d{\stackrel{\rightharpoonup}{r}}^{3}]}^{2}}},$$We should note that Eq. (2) explains the differences in contrast mechanisms for PAT and US. In PAT, the average absorption cross section is always greater than zero ($\overline{C}>0$), and the correlated power dominates the uncorrelated power for large $\rho $. As a consequence, the PA signal amplitude is approximately proportional to the optical absorption coefficient ${\mu}_{a}=\overline{C}\rho $ of the absorbing target. Similarly, when imaging a soft-tissue–bone interface with US from the soft-tissue side, the average backscattering coefficient is always positive ($\overline{C}>0$) and the correlated power dominates. When imaging soft tissue structures with US, however, researchers usually assume that the average backscattering coefficient^{17} is zero ($\overline{C}=0$), due to the way the scattered signals can be attributed to fluctuations in acoustic properties relative to the mean. Substituting $\overline{C}=0$ into Eq. (1) nullifies the second term on the right-hand side:

## Eq. (3)

$$\langle P(t)\rangle ={\sigma}_{\overline{C}}^{2}\rho {\int}_{V}{[\tilde{h}(\stackrel{\rightharpoonup}{r},t-|\stackrel{\rightharpoonup}{r}|/c)]}^{2}d{\stackrel{\rightharpoonup}{r}}^{3}\mathrm{.}$$As a result, US clearly relies on speckle contrast when imaging soft tissues. In the presence of multiple types of particles, Eq. (1) can be generalized to:

## Eq. (4)

$$\langle P(t)\rangle ={\int}_{V}{[\tilde{h}(\stackrel{\rightharpoonup}{r},t-|\stackrel{\rightharpoonup}{r}|/c)]}^{2}d{\stackrel{\rightharpoonup}{r}}^{3}[\sum _{i=1}^{M}({\overline{C}}_{i}^{2}+{\sigma}_{{C}_{i}}^{2}){\rho}_{i}]+{[{\int}_{V}\tilde{h}(\stackrel{\rightharpoonup}{r},t-|\stackrel{\rightharpoonup}{r}|/c)d{\stackrel{\rightharpoonup}{r}}^{3}]}^{2}{\left(\sum _{i=1}^{M}{\overline{C}}_{i}{\rho}_{i}\right)}^{2}\mathrm{.}$$We should clarify the definition of absorbers in PAT. In our model, unit absorbers must be 1) much smaller than the resolution cell, and 2) statistically independently and homogenously distributed in locations within a resolution cell. For example, when imaging blood vessels with a 5 MHz PAT system ($\sim 180\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ in resolution), we treat red blood cells (RBCs,$\sim 8\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ in diameter and $\sim 2\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ in thickness) as unit absorbers. Hemoglobin molecules, however, should not be treated as unit absorbers because their aggregation in RBCs differentiates their spatial distributions inside and outside RBCs and violates the condition in Eq. (2). When imaging water with megahertz PAT systems, individual water molecules can be defined as unit absorbers. We can also group every $K$ water molecules to form sparser super unit absorbers with a number density of ${\rho}_{\text{Sup}}=\rho /K$ as long as the dimension of the super unit absorber is much less than the resolution. The absorption cross section of the super unit absorber ${C}_{\text{Sup}}$ is the sum of the absorption cross sections ($C$) of the $K$ unit absorbers:

The location of the super unit absorber is the centroid of the locations of the unit absorbers. The super unit absorbers give rise to approximately the same image. Substituting Eq. (5) into Eq. (2), we have $({\overline{C}}_{\text{Sup}}^{2}+{\sigma}_{{C}_{\text{Sup}}}^{2}){\rho}_{\text{Sup}}=({\overline{C}}^{2}+{\sigma}_{C}^{2})\rho $ for the speckle term and ${\overline{C}}_{\text{Sup}}^{2}{\rho}_{\text{Sup}}^{2}={\overline{C}}^{2}{\rho}^{2}={\mu}_{a}^{2}$ for the boundary term. Therefore, the grouping leads to statistically equivalent results without violating our theory.

## 2.

## Method

## 2.1.

### Effects of Boundary Roughness on PAT Speckles: Simulation Studies

To analyze the effect of boundary roughness on PAT speckles, we simulate PAT of absorbing objects with boundaries having various degrees of roughness, which is quantified by the root-mean-squared (RMS) value ($\delta $) and the correlation length ($\xi $) of the boundary height.^{18} Our numerical phantoms are composed of a large number of absorbers, which are statistically independently and homogeneously distributed inside and outside the absorbing objects. The absorbers inside the absorbing objects have five times the average cross sections of those in the background. The boundary-profile functions of the absorbing objects are assumed to follow a stationary Gaussian stochastic process. The mean and the standard deviation of the stochastic process represent the mean boundary location and the standard deviation of the surface height (the RMS value $\delta $), respectively. The correlation function of the stochastic process is $G(r-{\stackrel{\rightharpoonup}{r}}^{\prime})=\mathrm{exp}(-|\stackrel{\rightharpoonup}{r}-{\stackrel{\rightharpoonup}{r}}^{\prime}{|}^{2}/{\xi}^{2})$, where $\stackrel{\rightharpoonup}{r}$ and ${\stackrel{\rightharpoonup}{r}}^{\prime}$ represent two position vectors on the mean boundary plane. Numerically, a boundary profile function is generated by convolving a zero-mean Gaussian distribution of random numbers for the surface height with the Gaussian correlation function.^{18}

The PAT simulation parameters are set to the same values as those of a custom-designed 512-element ring-array PAT system.^{19} Each ultrasonic transducer element is cylindrically focused in the elevational direction; therefore, an in-plane two-dimensional (2D) image can be reconstructed. The simulated mechanical-electrical impulse response (EIR) has a center frequency of 5 MHz (100% bandwidth). Using the Field II program,^{17}^{,}^{20} we calculate the spatiotemporal response of every ultrasonic transducer element due to all the absorbers in the imaging region. The PAT image is then reconstructed from these spatiotemporal responses.^{21} The in-plane spatial resolution, i.e., the full-width-at-half-maximum (FWHM) of the point spread function (PSF), is $\sim 180\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$.

## 2.2.

### PAT Speckle Statistics: Experimental Validations

We used the same ring-array based PAT system^{19} to experimentally study the statistics of PAT speckles. The illumination source was a tunable pulsed laser system based on an optical parametric oscillator (OPO; Vibrant HE 315I, Opotek, Inc.). The laser pulse had a repetition rate of 10 Hz and a pulse width of 5 ns. After collimation, the laser beam was homogenized through a ground glass before it reached the top surface of the tissue phantom. The absorbing objects were made of gelatin (10% gelatin by weight) mixed with graphite particles in various concentrations.

## 3.

## Results and Discussions

## 3.1.

### Effects of Boundary Roughness on PAT Speckles

The simulation results are tabulated by the RMS height $\delta $ and the correlation length $\xi $ in Fig. 1. As the RMS height $\delta $ increases, the boundary height fluctuates more; as the correlation length $\xi $ decreases, the position of each point on the boundary becomes less correlated with that of its neighbors; if $\delta =0$ or $\xi \to \infty $, the boundary becomes perfectly smooth [Fig. 1(a)]. Each cell of the table shows the true absorber distribution at the top and the corresponding reconstructed PAT image at the bottom. Segments of the boundaries of the absorbing objects within the dashed frames are extracted and shown as solid curves with a horizontal magnification of two times. The reconstructed PAT images show observable boundary buildups as well as interior speckles. In Fig. 1(a) ($\xi \to \infty $), and Fig. 1(b), 1(e), and 1(h) ($\xi =360\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$, which equals twice the in-plane resolution), the reconstructed and true boundaries agree well (correlation coefficient $\chi =0.92-0.95$) while the visibility of interior speckles ($V$) decreases slightly with increasing $\delta $. In Fig. 1(c), 1(f), and 1(i) ($\xi =180\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$, which equals the in-plane resolution), some of the features of the true boundaries cannot be recovered by the reconstructed boundaries ($\chi =0.73-0.77$). In Fig. 1(d), 1(g), and 1(j) ($\xi =30\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$, which equals $1/6$ of the in-plane resolution), the height fluctuations of the reconstructed and true boundaries become less correlated ($\chi =0.26-0.27$). Moreover, the amplitude of the reconstructed boundary vanishes with increasing $\delta $, increasing the visibility of interior speckles. In all the cases, however, the average powers of the interior speckles without normalization to the boundary signals are equal. Therefore, the variations in speckle visibilities are due to the changes in the amplitudes of the reconstructed boundaries.

In Fig. 2(a) and 2(b), the correlation coefficient ($\chi $) between the reconstructed and the true boundary profiles is quantified as functions of the correlation length $\xi $and the RMS height $\delta $. We found that $\chi $ depends solely on $\xi $ because the correlation coefficient between $\chi $ and $\delta $ is 0.08. If $\xi >180\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ (the in-plane resolution), all the boundary features can be resolved, and thus the reconstructed boundaries agree with the true boundaries. Conversely, if $\xi <180\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ (the in-plane resolution), features are too fine to be resolved. Therefore, even when the boundary fluctuations are small ($\delta =30\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$), the reconstructed boundaries do not agree with the true boundaries. In Fig. 2(c) and 2(d), the speckle visibility ($V$) is plotted as functions of the correlation length $\xi $ and the RMS height $\delta $. Because the interior speckle amplitudes follow Gaussian distribution, we define boundary features as the reconstructed image regions that have magnitudes more than three time of the standard deviation of the interior speckle amplitudes. We found that $V$ depends on both $\xi $ and $\delta $. In Fig. 2(c), $V$ decreases as $\xi $ increases when $\xi <180\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ (the in-plane resolution). If the correlation lengths are smaller than the in-plane resolution, shorter correlation lengths usually introduce more randomized phase-delay variations, especially when the RMS heights are larger than the in-plane resolution so that effective roughness presents. An extreme example is that if $\xi \ll 180\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$and $\delta \gg 180\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$ the reconstructed boundaries become fully developed speckles and the speckle visibility $V\to 1$. As $\xi $ increases further, the speckle visibilities $V$ at various $\delta $ gradually converge to the value of $V$ for the smooth boundary, because $\xi \to \infty $ indicates smooth boundaries regardless of $\delta $. In Fig. 2(d), if $\xi <180\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$, $V$ increases as $\delta $ increases, which introduces more randomized phase-delay variations. If $\xi \ge 180\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$, however, some reconstructed boundary segments [e.g., the paired arrows in Fig. 1(f)] have higher signal strengths than the reconstructed smooth flat boundaries. The geometric shapes of these features match the ring-shaped ultrasonic detection aperture, and therefore the phase variations of the PA partial waves from the surface absorbers are smaller than those from the smooth flat boundaries. As a result, the average boundary strengths of these boundaries are slightly higher than those of smooth boundaries, and thus the speckle visibilities are slightly lower.

The analysis above indicates that PAT speckles may appear in rare cases. For example, we may observe speckles when imaging melanoma, if the melanoma has rough boundaries whose correlation length is much smaller than and the RMS value is much greater than the imaging resolution. In contrast, when imaging blood vessels with most of the PAT systems (center frequencies ranging from 1 to 50 MHz), we observe only the boundary signals. Here, the RBCs are treated as individual absorbers. The RBC density in the blood of an average adult is $\sim 5\times {10}^{6}/{\mathrm{mm}}^{3}$. The blood vessel walls are considered smooth on the scale of the acoustic wavelengths. According to our analysis, the interior speckles are suppressed by the prominent boundary signals.

## 3.2.

### PAT Speckle Statistics

Photographs of the phantoms with high, medium, and low graphite particle concentrations and with smooth boundaries are shown in Fig. 3(a)–3(c). The corresponding PAT images are shown in Fig. 3(e)–3(g). As predicted in our previous study, we observed strong boundary buildups, which suppress the internal speckle patterns [Fig. 3(i)]. Also, the boundary features are most prominent at the highest particle concentration [Fig. 3(e)], and interior speckles become noticeable at the medium particle concentration [Fig. 3(f)] and become apparent at the lowest particle concentration [Fig. 3(g)]. The phantom with low particle concentration was imaged 10 times and 100 times, and the averaged reconstructed interior textures are shown in Fig. 3(j) and 3(k). Because both the phantom and the imaging system are stationary in each experiment, such averaging does not diminish the speckles but does reduce random noises. The similarity between the two interior images, with a correlation coefficient of 0.996, confirms that the interior texture is due to speckles rather than random noises. For comparison, a phantom with medium particle concentration and rough boundaries [Fig. 3(d)] was studied, where both the correlation length $\xi $ and the RMS height $\delta $ of the boundary profile are $\sim 60\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \mathrm{m}$. In the PAT images of the phantoms, the rough boundaries [Fig. 3(h)] produced 2.8 times weaker boundary amplitudes than the smooth boundaries at the same particle concentration [Fig. 3(f)].

We quantified the first order statistic of the PA speckles by plotting the histogram of the interior speckle amplitude (without the envelope detection) [Fig. 4(a)]. Since the coherent interference of ultrasonic waves can be described as a random walk process, the amplitude of the PA speckles follows a Gaussian distribution. The mean of the speckle amplitude is zero, because the PSF does not contain DC. The standard deviation of the Gaussian distribution ($\sigma $) which represents the square root of the average speckle power, is proportional to the product of the average particle absorption cross section $\overline{C}$ and the square root of the particle concentration ($\sqrt{\rho}$). In Fig. 4(b), we show that $\sigma $ is proportional to $\sqrt{\rho}$, while the boundary magnitude is proportional to $\rho $. As a consequence, the speckle visibility ($V$) is inversely proportional to the square root of the particle density $\sqrt{\rho}$ (15) for smooth boundary targets [Fig. 4(c)]. The second-order statistic of the speckles is shown in Fig. 4(d). In the classic speckle theory, the autocorrelation of the fully developed speckle pattern carries only the information of the system PSF rather than that of the target texture. The autocorrelation of the system PSF is shown in Fig. 4(e), which agrees with the autocorrelation of the speckle patterns [Fig. 4(f)].

## 4.

## Conclusions

In summary, we have studied the effect of boundary roughness on PAT speckles. The correlation coefficient ($\chi $) and the speckle visibility ($V$) were quantified as functions of the RMS value ($\delta $) and the correlation length ($\xi $) of the boundary height. If $\delta $ is much greater than and $\xi $ is much smaller than the imaging resolution, the reconstructed boundaries become fully developed speckles. In other words, speckle formation requires large uncorrelated height fluctuation within the resolution cell. We also experimentally studied the first- and second-order statistics of PAT speckles. The amplitude of the speckles follows a Gaussian distribution. The autocorrelation of the speckle patterns tracks that of the system point spread function. Although our analysis here was based on PA computed tomography, the linearity of PAT ensures that the conclusions hold for all PAT implementations.

## Acknowledgments

The authors would like to thank Li Li, Yan Liu, and Alejandro Garcia-Uribe for their helpful discussions. This work has been supported in part by National Institutes of Health Grant Nos. R01 EB000712, R01 EB008085, R01 CA134539, R01 CA157277, R01 EB010049, and U54 CA136398. L.W. has a financial interest in Microphotoacoustics, Inc. and Endra, Inc., which, however, did not support this work.