Retinal pigment epithelium (RPE) is a monolayer of pigmented cells located between photoreceptor outer segments and Bruch’s membrane. RPE has the function of transportation of nutrients and metabolic products between photoreceptors and blood^[1,2]. Thus, it plays an important role in the maintenance of retinal homeostasis and is implicated in a majority of retinal diseases due to the failure of RPE functions^[3–5]. An imaging method to examine RPE cells for retinal pathological analysis is needed.

Photoacoustic microscopy (PAM) is a noninvasive and label-free three-dimensional imaging modality based on the optical absorption property of biological tissues^[6–9], which has proved its capacity for in vivo imaging of blood vessels and RPE in fundus^[10–13]. To achieve the cellular or subcellular imaging level, PAM systems have been developed to achieve micrometer or submicrometer resolution and used for imaging red blood cells and pigment cells^[7,14,15]. Since the photoacoustic (PA) signals of RPE are mostly generated from the strong light absorption of the pigment components, especially the melanin, the PAM can display the cell morphology, and more importantly, indicate the melanin content in cells^[16], which is highly correlated with aging and retinal diseases such as age-related macular degeneration (AMD)^[5,17].

To acquire the RPE cell features such as morphology and component content, accurate cell detection and segmentation in cellular images are important. In previous work, several semi-automatic and automatic cell segmentation methods have been used for cell evaluation, and they acquired remarkable results in the images of corneal endothelium, photoreceptor, and RPE cells^[18–24], using two-photon fluorescence microscopy, confocal microscopy, and adaptive optics retinal imaging. As these imaging results cannot provide melanin content information of RPE cells, an automatic RPE cell segmentation algorithm for PAM imaging is necessary.

In our previous work^[15], RPE cells were imaged using dual modal PAM and optical coherence tomography (OCT). An ordinary region growing method was applied to segment RPE cells in PAM images by utilizing the connectivity of the relatively similar signal intensities within each cell. But the ordinary region growing method is limited as static predefined thresholds are required to avoid overgrowing of the segmented regions when applied directly to images^[25,26]. Since the signal intensities are not consistent between different cells in PAM images, no single static threshold is appropriate for the region growing process of each seed point. Thus, the previously adopted algorithm is unable to automatically segment all cells in PAM images.

In this Letter, to further realize automatic cell segmentation for PAM imaging and quantitatively studying cell features, we developed an improved region growing method by using the morphological characteristic of the RPE cell as a guide to dynamically update the growing threshold, so that the algorithm was capable of automatically traversing and identifying RPE cells in PAM images. Our automatic segmentation algorithm can potentially provide a quantitative analytical method for cell pathological research and disease detection at the cellular level using PAM imaging.

The experimental system for RPE cell imaging is a homemade submicrometer resolution PAM system, which is shown in Fig. 1. Details about the system description can be seen in our previous work^[15]. The lateral resolution of the PAM system is measured as 0.56 μm, and the irradiation energy on the samples is 10 nJ per pulse.

The region growing method for segmentation is based on an iterative approach to detect intracellular regions by adding connected points from the selected seed points in the image. Therefore, the accuracy of the cell segmentation relies on the positions of the seed points and growing thresholds to determine whether the neighborhood should be added to the region. In this Letter, we developed an algorithm by dynamically distributing seed points and adjusting segmentation thresholds to realize automatic RPE cell segmentation. The algorithm consists of three major steps, as shown in Fig. 2.

Step 1 is the seed point distribution. The seed points for region growing are initially uniformly spaced in the image. Some of the initial seed points are probably located in the intercellular space. These seed points will cause incorrect segmented cell regions, thus, the seed points need to be redistributed into intracellular regions. The seed point redistribution is implemented by first selecting a square region centered at the current seed position and then picking the point with the largest signal amplitude within the selected region as the new seed point. Since the PA signals are mostly induced by the strong light absorption of melanin within the RPE cells, the point with the largest signal amplitude corresponds to a point within the RPE cells if the neighborhood region is large enough to contain parts of the intracellular region.

In Fig. 3, the process of seed point distribution is illustrated in a PAM image of RPE. The initial distribution of seed points is shown as the red points in Fig. 3(a), and several points are located within the intercellular region. The dashed box in Fig. 3(a) represents the defined square region for the seed point in the center of the box. The size of the square region is 9×9pixels for the redistribution of each seed point. The blue point is the maximum value within the square region, which is also located in one RPE cell. The seed point is then modified to the blue point. After redistributing all of the initial seed points in Fig. 3(a), the redistributed seed points are shown as the blue points in Fig. 3(b), and the seed points are all in the intracellular regions for the following region growing. The density of the seed points and the range of neighborhood region can be modified based on the balance of segmentation accuracy and implementation time.

Step 2 is automatic segmentation. The segmentation is realized based on the region growing method using eight connected neighborhoods. First, an identical growing threshold Tg is used for region growing of each seed point and provides the segmented result as a binary image with one for the segmented region and zero for the rest. However, the obtained segmented region from this simple region growing will probably contain more than one cell due to overgrowing under an improper threshold. Thus, the minimum bounding rectangle (MBR) of the segmented region is calculated to judge the segmentation result and adjust the growing threshold Tg. Since the RPE cell has the morphology of a regular polygon, the MBR of the correctly segmented cell will have an aspect ratio (AR) close to one. When the segmented region overgrows to two or more cells, the AR will be much larger than one. Furthermore, if the segmented region consists of more than one cell, the region will appear as an irregular shape, which can cause significant difference between the area of the MBR and the segmented region. Thus, the judging condition can be written as AR=LMSM,R=ASAM,(1)S=AR<T1∧R<T2,(2)where LM and SM are the lengths of the long and short sides of the MBR, respectively, AR is the aspect ratio of the MBR, which represents the morphological characteristic of the segmented region, AS and AM are the areas of the segmented region and its MBR, respectively, R is the area ratio of the segmented region and its MBR, and T1 and T2 are two predefined judging thresholds. For positive judging result S, the cell segmentation for the current seed point is completed, and the next seed point will be loaded for segmentation. Otherwise, the growing threshold Tg will be decreased by a constant step St and used for a new segmentation iteration until the judging condition is satisfied.

Figure 4 illustrates the automatic segmentation process with the improved region growing method. Figures 4(a) and 4(e) are two PAM images of RPE cells with red arrows pointing to the selected cells for segmentation. The predefined judging thresholds T1 and T2 are chosen as 1.5 and 1.2 for each segmented region. Figures 4(b) and 4(f) are the intermediate results of the segmentation method with the MBR shown as the red boxes. Both intermediate results need to be re-segmented. The MBR in Fig. 4(b) shows a much larger AR than one when more than one cell is segmented. In Fig. 4(f), the AR does not have significant differences compared to the segmented result of only one cell, and the segmented region is not much larger than the calculated largest cell. In this case, the area ratio R in Eq. (1) between the segmented region and its MBR is important for further identification of the inaccurately segmented region. Therefore, the combination of judging conditions with AR and R enhances the accuracy for single cell identification and segmentation. The segmented results finally converge to the regions in Figs. 4(c) and 4(g), which match the original cells very well, as shown in Figs. 4(d) and 4(h).

Step 3 is post-processing. After segmentation for all seed points in the image is completed, a two-step post-processing is implemented. First, the segmented regions that are connected with the boundary of the image are removed. Secondly, since more than one seed point can be redistributed into an identical cell, repeated segmented results should also be removed. The cell repetition index is utilized for this purpose and defined as CIci/m=1−∑x,y|Ici(x,y)−Im(x,y)|Ici(x,y),(3)where Ici is the ith confirmed cell segmented binary image with the cell region of one and the rest zero, Im is the mth non-confirmed segmented binary image, and x, y are the pixel positions in the images. Cell repetition index CIci/m represents the difference between the mth result and the ith confirmed result by calculating the sum of the value difference between each corresponding pixel of the two images. The segmented results of different cells are rarely overlapped, and small CIci/m values can be obtained, while segmented regions of the same cell are largely overlapped, and the CIci/m values can be close to one. Thus, threshold T3 is used to judge whether the segmented regions are repeated. If the calculated cell repetition indexes all satisfy the condition of smaller than T3 after comparing it with the confirmed results, the mth result will be added into the confirmed cell sequence; otherwise, it will be removed as a repeated segmented result. Threshold T3 is chosen as 0.6 to remove repeated segmented results.

To validate the algorithm, a phantom study is designed to test the accuracy of the algorithm. The initial phantom image is shown in Fig. 5(a), which contains RPE-like structures created by generating a quasi-symmetric hexagonal array using MATLAB. The pixel value is one within the cell region and zero within the intercellular region. The dimensions of the image array and pixel sampling were similar to those in real PAM images.

Phantom images of various signal-to-noise ratios (SNRs) are utilized to test the algorithm. Since the PAM image is shown as the maximum-amplitude-projection (MAP) image in which each pixel corresponds to the maximum value of the PA signal at each imaging position, we first simulated the PA signal at each pixel of the initial phantom image, and then we added white Gaussian noises to form a noisy PA signal. The noisy phantom image can be acquired by getting the maximum value of the noisy simulated PA signal at each pixel position. The SNR of the noisy phantom image is defined as the ratio between the mean pixel value of the cell region and the intercellular background. A typical acquired RPE cell amplitude profile in a PAM image is shown in the top of Fig. 5(b), which is similar to the simulated cell amplitude profile with a 12 dB SNR. Simulated cell amplitude profiles under different SNRs are also illustrated in the bottom of Fig. 5(b), and the identification of adjacent cells is increasingly difficult as the SNR drops.

The segmentation algorithm was tested using phantom images under different SNRs with the normalized 5 and 3 dB phantom images shown in Figs. 5(c) and 5(e). The segmented result of the 5 dB phantom image shows clear and complete hexagonal cell morphology, while irregular and incomplete segmented cell regions appear when the SNR is reduced to 3 dB, as shown in Figs. 5(d) and 5(f). To quantify the algorithm performance, we made a statistic analysis of the segmented cell number and area under different noise levels. The cell number and area were acquired by calculating the number of segmented regions and the pixel number within each segmented region, respectively. The segmentation under each noise level was repeated four times to get the statistical result. As shown in Figs. 5(g) and 5(h), when the SNR of the phantom image is 5 dB or larger, the cell number and area remain consistent and show no significant differences. When the SNR is 4 dB, the segmented cell number is still constant, while cell area becomes smaller, which means incomplete cell regions are segmented due to the increased noise level. When the SNR is below 4 dB, the segmented cell number starts to change, and the cell area becomes even smaller and more inconsistent, as the algorithm may miss several cells or segment more than one region within one single cell in phantom images with a very low SNR. Since the SNR is similar to 12 dB in real PAM image, as shown in Fig. 5(b), the automatic segmentation algorithm has the capability of correctly segmenting most RPE cells and obtaining reliable cell information.

A piece of the RPE layer stripped from a normal porcine eye was used for PAM imaging. The PAM image and the subsequent segmentation and quantitative results of RPE cells are shown in Fig. 6. The PAM image is the MAP image in which each pixel value is calculated as the maximum value of the acquired PA signal at the corresponding scanning position. The segmentation result of all of the complete cells in Fig. 6(a) is shown in Fig. 6(b). Most segmented cells show typical morphological characteristics of the RPE cell and match the original cell regions, but there are still some cells that may not be segmented completely due to the low SNR within a local region of the PAM image.

After segmentation, the cell area is acquired by calculating the area of the segmented region. The cell intensity that mainly reflects the melanin content in the RPE cell is calculated by the sum of the PA signal amplitudes within the segmented cell region. In Fig. 6(c), a statistical analysis of the cell area is presented. The area distribution is centered at 200μm2, which is similar to the result of the normal human eye^[23]. The relationship between the cell area and intensity is shown in Fig. 6(d). Two regions of the stripped RPE were imaged and segmented for quantitative analysis under the same experimental conditions and segmentation parameters. The cell area and intensity are linearly related for normal porcine RPE cells, shown as the fitted red and blue lines (R2>0.7). The two fitted lines show similar slopes with a mutual difference of less than 5%, which proves the consistency of the PAM imaging and segmentation algorithm. Since the RPE cell suffers from morphology changes and melanosomes leaking at the early state of the diseases such as AMD^[5], the calculated distributions of cell area, intensity, and their correlation can be promising indicators for detecting cell abnormality caused by diseases.

In conclusion, we develop an improved region growing method for RPE cell analysis in PAM images to realize completely automatic cell segmentation and feature calculation of both cell area and component content, which can have the potential of detecting cell abnormality and early retinal diseases in a future study. Furthermore, this method may also be applied in other cytological studies when the judging condition in the algorithm is adjusted correspondingly.