Automatic threshold selection algorithm to distinguish a tissue chromophore from the background in photoacoustic imaging

: The adaptive matched filter (AMF) is a method widely used in spectral unmixing to classify different tissue chromophores in photoacoustic images. However, a threshold needs to be applied to the AMF detection image to distinguish the desired tissue chromophores from the background. In this study, we propose an automatic threshold selection (ATS) algorithm capable of differentiating a target from the background, based on the features of the AMF detection image. The mean difference between the estimated thickness, using the ATS algorithm, and the known values was 0.17 SD (0.24) mm for the phantom inclusions and -0.05 SD (0.21) mm for the tissue samples of malignant melanoma. The evaluation shows that the thickness and the width of the phantom inclusions and the tumors can be estimated using AMF in an automatic way after applying the ATS algorithm.


Introduction
Multispectral photoacoustic (PA) imaging is a rapidly growing modality that can provide highresolution target detection based on the optical characteristics of tissue [1][2][3]. PA imaging uses ultrasound beamforming to locate the source of the acoustic waves generated by the thermoelastic response following the absorption of nanosecond laser pulses in biological tissue. PA imaging provides more specific information on individual tissue components than ultrasound imaging, together with high spatial resolution. PA imaging is less sensitive to scattering than optical imaging, and is capable of imaging endogenous tissue chromophores deeper in the tissue [4][5][6]. The optical absorption of photons, and thus the PA response, will change depending on the wavelength of the transmitted laser pulse. This provides a specific spectral signature for each tissue chromophore [7,8]. However, this also results in wavelength-dependent light attenuation leading to changes in both the amplitude and morphology of the PA response, called "spectral coloring" [9,10], making complex spectral processing necessary.
Spectral unmixing is a unique approach capable of classifying various optical absorbers based on their spectral signatures. Spectral unmixing has shown high target detection ability in both extrinsic and intrinsic optical absorbers [11][12][13][14][15]. It has been used in the visualization of the vascular system [16,17], for blood oxygen saturation mapping [18,19], the detection of human carotid plaques [20,21], and the detection of basal cell carcinomas [22,23]. Various spectral unmixing strategies have been developed to distinguish chromophores of interest from the absorbing background, such as linear spectral fitting [24][25][26], statistical sub-pixel detection [27,28], and blind source unmixing methods [16,29]. The performance of spectral unmixing methods depends on the excitation wavelength range, the tissue chromophore spectra, the number of chromophores, and the ability to separate intrinsic and extrinsic chromophores from the background tissue [30,31].
The adaptive matched filter (AMF) [32,33] is a conventional statistical sub-pixel detection method used for spectral unmixing, and has shown encouraging target differentiation results in both phantoms and clinical studies. This method has shown promising results in hyperspectral imaging when the target spectral signature is known, but precise information on the background optical absorbers is lacking [34,35]. It has been shown in a previous study that AMF target detection outperforms common linear unmixing algorithms in the presence of sparse target distribution in the background [36]. AMF estimates the probability of each pixel being in the target area based on the statistical characteristics of the background. A Gaussian distribution is used to model the background, as this can reduce the effect of spectral coloring [37]. Applying AMF to PA images at different wavelengths results in a detection image, in which the amplitude is related to target probability. However, exact discrimination between target and background pixels is challenging.
Histological examination in clinical studies has shown that greater thickness of the skin tumor is correlated with a greater risk of metastases [38]. In many cases, precise information is required on the target distribution, the target thickness in the background, and the number of target pixels (area, volume). In order to differentiate the target from the background in PA imaging and delineate its borders, a threshold needs to be applied to the AMF detection outcome. This threshold can be set by manual fine-tuning by an operator. However, this requires highly-skilled operators, and the results may differ between operators. Another approach is to define a fixed threshold for all the images, but this may lead to false detection in some of the outcomes depending on the PA measurement conditions.
In this study, we address the complexity of threshold selection by developing an automatic threshold selection (ATS) algorithm that can be applied to AMF detection images. The performance of the proposed technique was investigated on a tissue-mimicking phantom with inclusions of different sizes. The feasibility of its use in tissue was also studied on ex vivo human tissue samples of malignant melanoma, including comparisons of the thickness and width of the tumors determined histologically.

Materials and methods
The chromophore distinguish technique, included three steps: 1) preprocessing of multispectral PA images, 2) application of the AMF technique, and 3) the implementation of an ATS algorithm. The processing steps were carried out using MATLAB (The MathWorks, Inc., Natick, MA, USA).

Preprocessing
Preprocessing was carried out as follows. Pixels with an amplitude lower than 2% of the maximum amplitude were discarded. The pixel spectra were smoothed in the wavelength direction using a moving average filter with a window size of 3 to suppress noise, while capturing the important patterns in the data.

Application of the adaptive matched filter
AMF computes a pixel-wise detection value to determine the similarity between the spectrum in each pixel and the target spectrum. A higher output value corresponds to a higher similarity to the target spectrum [39]. AMF is defined by where D is the AMF detection value for each pixel at position r in the image, µ t is the target reference spectrum, µ b is the background mean spectrum, Σ is the background covariance matrix, T is the matrix transpose operation, and x r corresponds to the multispectral data at position r. D is a measure of the probability that a spectrum in a pixel belongs to the target. By applying a threshold to D, pixels can be divided into either target or background. In this study, AMF resulted in an m × n image with different detection values, where m and n are the numbers of pixels in the vertical and horizontal directions, respectively. The target reference spectrum, µ t , was defined as the mean of the target spectra from the multiwavelength PA images. In cases where the target was sparsely present, all image pixels were selected to estimate the background covariance matrix Σ and the background mean spectrum µ b , otherwise, the background spectra were extracted from the pixels known to be background.

The algorithm for automatic threshold selection
In order to quantify the number of pixels within the target area precisely, a threshold needs to be set. In this work, an algorithm for automatic threshold selection was applied to an AMF detection image to separate the target pixels from the background automatically. The algorithm includes several steps, as outlined in Fig. 1. In this approach, decreasing a threshold from the maximum detection value found by AMF, denoted D max , to the minimum detection value, denoted D min , gives S different images, in which the pixel values lower and higher than the thresholds are set to 0 and 1, respectively. S is the number of steps used in the algorithm. The step size s is obtained by In this study, S was set to 100. When the threshold is decreasing, an increasing number of pixels are detected as target pixels. During this process, several target pixels are connected and formed to a group, in this study denoted as a connected component. This happens at several locations in the image, and as the threshold is decreased, the ATS algorithm calculates two features in each step: 1) the number of connected components, and 2) the number of pixels in the largest connected component. We found that these features, as the threshold decreases stepwise, results in a behaviour that can be used to automatically choose a suitable threshold. The MATLAB functions "bwconncomp" and "regionprops" were used to obtain the connected components. Once all the steps have been completed, this result in two feature vectors: the number of connected components for each threshold, denoted G = [G 1 , G 2 , · · · G L ], and the number of pixels in the largest connected component for each threshold, denoted H = [H 1 , H 2 , · · · H L ]. Different shapes of the curve describing the number of connected components, while decreasing the threshold, were observed. In order to obtain robust measurement, these curves needed to be smoothed. We found that a two-term Gaussian model (F) fitted to G gave a simple constrained shape for automatic measurement, while still yielding a good fit. Figure 2 shows four examples of the G and F obtained from the tissue-mimicking phantom ( Fig. 2(a,b)), and the ex vivo human tissue samples of malignant melanoma (Fig. 2(c,d)). Figure 3 shows F and H while increasing the number of steps. A distinct antegrade can be seen in F, due to an increasing number of sparsely connected components, which is followed by a distinct retrograde, resulting from the merging of various sparsely connected components. As the number of steps increases, the number of connected components starts to increase (point P1), meaning that the target has started to expand. The maximum number of connected components, is denoted P2, see Fig. 3(a). However, at P2, some of the background pixels have erroneously merged with the target. Hence, a potential stopping point can be found between P1 and P2, denoted P3. P3 is determined as the point of the highest abrupt change in the range [P1, P2], as a very high abrupt change in this range suggests that the number of connected components increases sharply due to false background pixels detection ( Fig. 3(b)). The MATLAB function "findchangepts" is used to find the highest abrupt change. This function returns the index at which the mean and slope of the data changes most significantly. Point P3 thus provides an estimate of the threshold. In order to determine a more accurate stopping point, H was used to follow the changes in size of the largest connected component (Fig. 3(c)). If an abrupt change is found when taking the adjacent step numbers to P3 in H, then P3 is replaced by T.
Step numbers higher than T cause background pixels to be falsely merged with the targets. In this study, a range of 10 pixels adjacent to P3 was selected to search for T (Fig. 3(d)). If a sudden change in H is not observed, then P3 is the selected threshold, and is denoted T. Pixels in the AMF detection image that exceed the threshold corresponding to step T are classified as the target.

Photoacoustic imaging system
PA imaging was performed using a Vevo LAZR-X imaging platform (FUJIFILM VisualSonics Inc., Toronto, ON, Canada). Two planar laser beams at the sides of the ultrasound transducer were used to illuminate the surface of the tissue-mimicking phantom and ex vivo tissue samples, and the PA data were acquired using an MX250 ultrasound linear array transducer. Multispectral data were recorded in the wavelength range from 680 nm to 970 nm, in steps of 5 nm.

Tissue-mimicking phantom
Target detection was first performed in a tissue-mimicking phantom in order to test the ability of the algorithm to detect targets of various sizes. The phantom was made of styrene-ethylene/butylenestyrene [40], with the dimensions 30 × 20 × 5 mm, and included four different homogeneous inclusions as the targets, as shown in Fig. 4(a). The inclusions were hemispherical with radii of 3 mm, 2 mm, 1 mm, and 1/2 mm. A green pigment was used as the background, and a mixture of green and black pigments was used in the inclusions. Two separate molds were used to make the background and the inclusions. They were then assembled in the third mold, and placed in an oven, with the temperature set to 130°C, to fuse the parts together. The phantom was scanned in 3 dimensions using a stepping motor with a step size of 0.5 mm ( Fig. 4(b,c)). The photoacoustic spectral signatures of the targets and the background are shown in Fig. 4(d).

Ex vivo tissue samples
Ex vivo multiwavelength PA imaging was carried out on seven maligmant melanoma samples from human skin to investigate the performance of the ATS algorithm on tissue. 3-D PA scanning was carried out as described above for the phantom. The mean spectrum of all parts confirmed to be tumor was used as the target signature. The PA images were analyzed using the AMF spectral unmixing method, and the results used as input for the ATS algorithm to distinguish tumor pixels from the healthy tissue without the need for a highly skilled operator. A comparison between ATS and a number of fixed thresholds was also conducted. These tumors have previously been published in a clinical study on PA imaging and melanoma. Figure 5 shows the detection images obtained from the four inclusions in the phantom, with and without the ATS algorithm. A clear border can be seen between the inclusions and the background in the images when ATS is applied, while it is lacking when ATS was not applied. Despite the fact that the inclusions were homogeneous, the deeper parts of the inclusions do not show higher AMF values than the surface of the layer. This is due to the spectral coloring and the high absorption of the inclusion (Fig. 5(a)). Figure 5(b) shows the ATS result overlaid on the ultrasound image, where the thickness of the inclusions estimated using ATS, and approximate extent of the inclusions are indicated. Despite the heterogeneous AMF detection images without ATS, border of the target is clear in the images with ATS. The AMF detection image size (pixel height × pixel width) used in the phantom study was 100 × 400, and the average ATS algorithm run time was 1.2 seconds (Processor: Intel Core i7 CPU @ 3.40 GHz). A scatter plot of the thickness estimated when applying ATS versus the known thickness of the four inclusions is shown in Fig. 6(a). A root mean squared error of 0.26 mm was obtained. Figure 6(b) shows a Bland-Altman plot for the four different inclusions. No data were outside the lower or upper confidence bounds, and the mean difference was 0.17 SD (0.24) mm.

Ex vivo tissue samples results
The results regarding skin tumor detection using AMF with and without ATS are shown in Fig. 7. It can be seen from Fig. 7(a), which shows the results of AMF without ATS, that the border of the tumor cannot be distinguished from the surrounding healthy tissue, and that estimates of the tumor thickness would not be accurate with a fixed threshold for all samples, especially in Subjects 2, 3, and 4. the results obtained when applying ATS are shown in Fig. 7(b), where the part of the skin detected as tumor is overlaid on the ultrasound images. The greatest tumor thickness is indicated in the image for each subject. Figure 7(c) shows the results of histological examination, including the manually estimated tumor thickness and dimensions. The AMF detection image size (pixel height × pixel width) used in the ex vivo study was 656 × 512, and the average ATS algorithm run time was 2.3 seconds (Processor: Intel Core i7 CPU @ 3.40 GHz).
A scatter plot and a Bland-Altman plot for the thickness of the tumors obtained from ex vivo tissue samples are shown in Fig. 8. The root mean squared error was 0.19 mm, and the mean difference was -0.05 SD (0.21) mm. As for the phantom, no estimate was outside of the lower or upper confidence bounds. Figure 9 compares the width of the tumors estimated using ATS with the results obtained from histological examination. One of the samples was outside of the upper confidence bound (Subject 6). The mean difference between the width estimated with ATS and the histologically determined width was 0.50 SD (3.52) mm.
The statistical results for the thickness and width estimates for each ex vivo tissue sample of malignant melanoma are summarized in Table 1. The statistical results concerning the comparison of the thickness measurements using ATS in both the phantom and ex vivo tissue samples, to the known values and histological results are given in Table 2. Figure 10 compares the absolute difference of the tumor thickness and width estimated using ATS and a number of fixed thresholds (0.5, 0.4, 0.3, and 0.2) in ex vivo human tissue samples of malignant melanoma. The median of the absolute error of tumor thickness and width measurements using the ATS were lower compared to using a fixed threshold.

Discussion
In PA imaging, a threshold needs to be applied to the AMF detection image in order to differentiate the target from the background and delineate the target area. The resulting delineation of target spectra can be an important criterion in tissue diagnosis, for instance when characterizing skin tumors [41]. To the best of our knowledge, this is the first time AMF detection image features were used to develop an automatic threshold selection algorithm that could be applied to estimate the thickness and extent of the target. The validity of ATS was demonstrated using multispectral PA imaging of a tissue-mimicking phantom and ex vivo tissue samples of malignant melanomas. An AMF spectral unmixing approach can be employed when the target spectrum is known but there is no reliable information on the background optical absorbers [34], as was the case in the current study. In a previous study, it has been shown that the AMF suppresses the effect of wavelength-dependent fluence variation by modelling the background and target spectral variability [36]. Therefore, in this first study on the ATS algorithm, wavelength-dependent fluence variation was not taken into account explicitly. However, future studies need to investigate the effect of this assumption, especially for deep tissue. The ATS algorithm could potentially be applied to spectral unmixing methods other than AMF, however, this requires further investigation.
Some difficulties were encountered when making the phantom. The background material and the inclusions were made separately and then assembled, and placed in an oven to allow them to fuse together. During this process, parts of the edge of some inclusions melted into the background, especially in inclusion 4, which was very small and difficult to make. In addition, we tried to analyze the deepest slice of each inclusion using 3-D scanning in order to compare the estimated thickness to the known value. However, the spatial resolution of the 3-D scanning was 0.5 mm, which is close to the thickness of the smallest inclusion. This caused problems in determining the deepest part of the smaller inclusions (3 and 4). In 3 of the 4 inclusions, the thickness estimated with ATS was close to the known thickness.
The purpose of evaluating ex vivo tissue samples was mainly to demonstrate the feasibility of applying ATS to human tissue. The alignment between the histological images and the images from the PA measurements was not exact, but attempts were made to analyze the PA images with the closest slice among the histological images. As no real ground truth was known, the histological results were used as measures of tumor thickness and width. In Subjects 3 and 4, the thickness estimated with ATS was greater than the thickness obtained from the histological images. This may be due to mismatch between the histological specimens and the analyzed slice. The median of the absolute difference of the estimated tumor thickness and width with ATS were lower than using a fixed threshold applied to the AMF detection images. Applying ATS in combination with AMF led to better agreement in tumor thickness and width measurement compared to thickness and width estimates without ATS. Since the ATS algorithm is applied to the AMF detection images, in samples where AMF was not able to include all parts of the tumor, the performance of the ATS impaired accordingly (Subject 6). The ATS algorithm is based on connectivity analysis, and in samples where the target area is not connected (Subject 6) further analysis is needed to estimate the width of the tumor. Also, in this study, there were no complex morphologies or boundaries, and the ATS performance on such samples needs to be evaluated in future studies.

Conclusion
We have presented an automatic threshold selection algorithm capable of differentiating the targets from the background in AMF detection images. The performance of the ATS was evaluated using a tissue-mimicking phantom. Excellent agreement was found between the estimated thickness of the inclusions and the known values with a mean difference of less than 0.2 mm. The feasibility of using the ATS algorithm to differentiate malignant melanoma from healthy tissue in ex vivo human tissue samples of skin tumors was investigated. The ATS showed better performance compared to using fixed threshold values. Further evaluation is needed, however, we envisage that the algorithm can be an important tool in future clinical investigations.