Machine Learning Analysis of Human Skin by Optoacoustic Mesoscopy for Automated Extraction of Psoriasis and Aging Biomarkers

Ultra-wideband raster-scan optoacoustic mesoscopy (RSOM) is a novel modality that has demonstrated unprecedented ability to visualize epidermal and dermal structures in-vivo. However, an automatic and quantitative analysis of three-dimensional RSOM datasets remains unexplored. In this work we present our framework: Deep Learning RSOM Analysis Pipeline (DeepRAP), to analyze and quantify morphological skin features recorded by RSOM and extract imaging biomarkers for disease characterization. DeepRAP uses a multi-network segmentation strategy based on convolutional neural networks with transfer learning. This strategy enabled the automatic recognition of skin layers and subsequent segmentation of dermal microvasculature with an accuracy equivalent to human assessment. DeepRAP was validated against manual segmentation on 25 psoriasis patients under treatment and our biomarker extraction was shown to characterize disease severity and progression well with a strong correlation to physician evaluation and histology. In a unique validation experiment, we applied DeepRAP in a time series sequence of occlusion-induced hyperemia from 10 healthy volunteers. We observe how the biomarkers decrease and recover during the occlusion and release process, demonstrating accurate performance and reproducibility of DeepRAP. Furthermore, we analyzed a cohort of 75 volunteers and defined a relationship between aging and microvascular features in-vivo. More precisely, this study revealed that fine microvascular features in the dermal layer have the strongest correlation to age. The ability of our newly developed framework to enable the rapid study of human skin morphology and microvasculature in-vivo promises to replace biopsy studies, increasing the translational potential of RSOM.


I. INTRODUCTION
N ON-INVASIVE and quantitative in-vivo assessment of skin features, including the microvasculature, carries significant potential for diagnostics and disease monitoring in a number of pathologies [1], [2].However, the use of non-invasive observations is currently limited by the tools available.Consequently, most of our knowledge on the implication of morphological and microvascular skin features in various skin and systemic diseases are based on histological analysis of biopsied skin samples ex-vivo [3], [4], [5], [6], [7].For example, psoriasis leads to epidermal thickening, capillary elongation and increased dermal vascularization [4], [5], [8].Aging, diabetes, and cardiovascular disease lead to changes in subcutaneous microvascular morphology and function [9], [10], [11], [12], [13].While histological sampling has shed light into these relationships, it is an invasive procedure associated with pain and risk of infection.In addition, biopsies are very laborious and costly, making them undesirable for routine examinations in comparison to diagnostic and theranostic tests that take into account microvascular alterations [14], [15].
While imaging techniques can be used to analyze skin non-invasively [1], [15], [16], many of the current methods do not offer the fine detail of histological analysis.For example, dermoscopy visualizes only the skin surface and is not appropriate for retrieving skin features under the epidermis due to the strong photon scattering of skin [15], [17].Tissue sectioning microscopy methods offer technology that can reduce the effects of photon scattering but still have limitations in skin imaging [1], [15], [16].Confocal microscopy requires high optical energy per volume element and only reaches a few tens of microns deep [18].Optical coherence tomography (OCT) [19], [20] and raster-scan optoacoustic mesoscopy (RSOM) reach deeper into the skin and can also extract vascular features.However, the optoacoustic method has deeper penetration and provides stronger contrast from the vasculature, holding the potential to become a widespread method for the study of skin morphology and function.It has been recently shown that RSOM is the only method available today that can offer three-dimensional skin images with virtually isotropic resolution and highly detailed crosssectional images [7], [8], [21].Moreover, the different chromophores that can be visualized at various wavelengths render RSOM superior over other methods in terms of the functional contrast visualized [21].
RSOM has been already employed in human studies to quantify psoriasis burden and remission due to treatment [3], [8], visualize vasculature associated with melanoma formation [7] or capture functional skin characteristics in response to heating [22].In addition, structural and functional imaging features derived from RSOM have been demonstrated to be objective disease markers to assess and stratify the severity of atopic dermatitis and were found to correlate well with conventional metrics [23], [24], [25].Depending on the implementation of RSOM, the method is able to acquire super broad bandwidth optoacoustic signals from ten up to more than a hundred MHz, resulting in resolutions down to tens of micrometers or better and penetrations of up to several millimeters [3], [7], [8], [21], [22], [24], [26], [27].Nevertheless, the laborious nature of processing RSOM images does not yet allow routine dissemination of the modality in clinical settings.Firstly, manual segmentation and simple intensity-threshold based approaches applied today are laborious, user-dependent and do not allow for testing of large patient cohorts in a time efficient manner [3], [8], [26].Moreover, user-dependent operation may introduce errors.Secondly, the fine structures of RSOM images correlated to high frequency optoacoustic signals typically attain much lower intensity compared to large skin structures reconstructed from low frequency signals, which may result in biases in image analytics especially when applying intensity thresholds for image segmentation [28].Thirdly, due to light attenuation or substantial variations in image contrast at different skin depths, there may be depth-dependent intensity loss in the image, which is an image feature that cannot be addressed with simple filtering operations.
To explore the full potential of RSOM and address these limitations, we employed deep learning as a tool to automate and improve the segmentation accuracy and to extract features and imaging biomarkers automatically.Convolutional neural networks (CNNs) represent the state-of-the-art for pixel classification (segmentation) in general computer vision and medical imaging modalities such as magnetic resonance imaging (MRI) and computed tomography (CT) [29], [30], [31].For vascular structures, the DeepVesselNet architecture [32], [33] and topology-preserving loss functions [34], [35], [36], [37] have been developed to automatically segment 3D vessel structures from MRI, CT and optoacoustic tomography images [38], achieving excellent segmentation performance.In this study, we demonstrate for the first time the We demonstrate DeepRAP by applying it to automatically segment RSOM images obtained from 25 psoriasis patients under treatment.The quantified RSOM biomarker is then applied to characterize the disease severity and progression, showing excellent agreement with manual segmentation and histology.In addition, we test DeepRAP using a more challenging problem, namely the assessment of cutaneous microvascular endothelial function by analyzing a sequence of RSOM volumetric images acquired during the post-occlusive reactive hyperemia process, i.e. an image sequence with significant contrast variations.Results show that DeepRAP accurately captures and quantifies the strong dynamic changes of skin microvasculature features, in higher detail and accuracy compared to Laser Doppler flowmetry or tissue spectrometry.We found that DeepRAP performs well, even at varying signal intensities due to tissue inhomogeneity at different skin depths or from different skin conditions.Having validated DeepRAP in datasets with known performance, we applied it to explore the rate of microvasculature change as a function of age in a group of 75 healthy volunteers.DeepRAP extracted five vessel features, which were examined for their relationship to age progression.The analysis indicates that small vessels in the 10-40 micrometers range were most prominently affected by age, with a reduction rate that appeared most prominent in the 20-65 years' age range.The combination of RSOM and DeepRAP analysis presents an attractive solution to image and quantify morphology and functional changes in the skin, with the potential to improve diagnostic and prognostic applications for skin and circulatory pathologies..

A. RSOM Imaging and Image Reconstruction
We employed an in-house RSOM imaging system, which was introduced in our previous work [3], [39].Illumination was provided by a pulsed laser at a wavelength of 532 nm.The repetition rate of the laser was 1 kHz, yielding an optical fluence of 0.375 mJ/cm 2 , which is far below the safety limit according to the American National Standards for Safe Use of Lasers in humans (20 mJ/cm 2 ) [40].Before each scan, the skin was cleaned with alcohol wipes.Both the patients and the operators used appropriate goggles for laser safety reasons.Each patient was scanned with an imaging field of view of 4 × 2 mm 2 , a step size of 7.5 µm in the fast axis (X axis), and a step-size of 15 µm in the slow axis (Y axis).Z is the depth direction.Each RSOM scan lasted approximately 70 s.We first applied motion correction algorithms to minimize motion-related artifacts in every RSOM scan before reconstruction [41].Then, acquired RSOM signals were divided into two frequency bands, 10-40 MHz (low) and 40-120 MHz (high), for the 10-120 MHz bandwidth.Signals in the two different bands were independently reconstructed.Reconstructions were based on beam-forming algorithms that generated three-dimensional images [39].The reconstruction algorithm was accelerated by parallel computing on a graphics processing unit (GPU) and improved by incorporating the spatial sensitivity field of the detector as a weighting matrix.The reconstruction time of one bandwidth took about 5 minutes with the voxel size of the reconstruction grid at 12 µm × 12 µm × 3 µm.The two reconstructed images R_low and R_high corresponded to the low-and high-frequency bands.A composite image was constructed by fusing R_low into the red channel and R_high into the green channel of an RGB image [3].The detailed process has been introduced in our previous work [3].The RSOM images were rendered by taking the MIPs of the reconstructed images along the slow axis or the depth direction.

B. Volunteer and Patient Studies
Twenty psoriasis patients with Psoriasis Area Severity Index (PASI) values from 1 to 7 were imaged following approval from the Ethics Committee of the Technical University of Munich.In addition, 5 psoriasis patients were measured during conventional inpatient treatment consisting of topical descaling, anti-inflammatory therapy (by means of salicylate Vaseline, topical corticosteroids and dithranol) and simultaneous phototherapy (311 nm Narrowband UVB (NB-UVB) or psoralen-UVA (PUVA)).RSOM scans were recorded at different time points on the same skin location.The detailed information of the study has been reported in our previous work [8].
For the hyperemia experiment, 10 healthy volunteers with a mean age of 33 were recruited following approval from the Ethics Committee of the Technical University of Munich.The RSOM scanning head was positioned at an area of about 5 cm to the wrist.A clinical use pneumatic cuff was placed at the level of the upper arm (i.e., distal to the site of brachial artery measurement) and controlled by an experienced operator.The hyperemia measurement took a total of 9 minutes, including: 2 minutes baseline, 4 minutes cuff on (cuff pressure was inflated to 220 mmHg), and 3 minutes cuff off (cuff deflation).To visualize the skin microvessels during the 9-minute cuff measurement, 9 RSOM 3D scans were recorded every minute in an area of 4 mm × 2 mm.
Healthy volunteers were recruited following approval from the Ethics Committee of the Technical University of Munich.All participants gave written informed consent before the planned RSOM examination.In total, 75 healthy volunteers, with ages ranging from 25-65 years, were scanned.The volunteers were divided into three age groups: I (n=24, 29.5 ± 3.5 years, in the range of 20 to 35 years), II (n=28, 41.1 ± 4.1 years in the range of 36 to 50 years) and III (n=23, 59.5 ± 4.9 years in the range of 51 to 65 years).Each individual was scanned at a region of interest (ROI, 4×2 mm) over the pretibial region of the distal lower limb.RSOM data quality was evaluated based on our previously developed RSOM quality evaluation approach and low-quality data was excluded [42].

C. Skin Layer Segmentation
Our network architecture is a U-Net with a depth of 5 and dropout in the second to fifth up-convolution, based on the work introduced by Gerl et al. [43].This architecture uses an encoder-decoder structure wherein the encoder uses blocks of convolution and max-pooling operations to encode the image information into a compressed feature space.The decoder takes these encodings from the feature space as input and uses up-convolutions and regular convolutions to decompress the information into our desired representation.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Finally, with residual connections, the output of each block during the encoding is fed forward to the corresponding decoder block.
For the layer segmentation model, 28 RSOM volumes from volunteers and 15 RSOM volumes from psoriasis patients are used for training, 7 RSOM volumes from volunteers and 3 RSOM volumes from psoriasis patients are used for validation and 8 RSOM volumes from volunteers and 2 RSOM volumes from psoriasis patients are used for testing.We split each volume into 2D slides from both sides, resulting in a training set of 8550 2D samples with a pixel resolution of 333 × 550.We used the following data augmentations: (1) random Z rescale (rescaling of the epidermis region by a random factor between 0.6 and 1.4); (2) random Z shift (shifting the whole volume in the z direction by a random value between −75 and 100); (3) random mirror; and (4) intensity transform (a piecewise linear intensity rescaling).We performed a five-fold cross validation.As a loss function, we used binary cross entropy.Furthermore, we used the Adam optimizer and trained our models on one P5000 GPU.We achieved a quantitative segmentation performance of 84.26 ± 8.22 in Dice score and 73.63 ± 11.71 in IoU.We then calculated two primitive skin layer features.First, we extracted the average epidermis width in µm.Here we divide the total epidermis volume by the fixed RSOM X and Y dimensions (333 and 171).
To enhance signal density in low-signal areas of the epidermis and to aid segmentation, we employed a 1D sliding window maximum filter along the direction perpendicular to the 2D slice fed into the U-Net.We explored this post-processing step in an experiment, where we applied different filter lengths to a set of 10 RSOM volumes.We then calculated the resulting segmentation performance in segmentation metrics and identified disconnected components in the samples.For a set of 10 RSOM volumes, the segmentation metrics (i.e., Dice, Precision, Recall) for different filter lengths were comparable as shown in Table I.However, at a length of five, we saw a minimum in the number of disconnected components of the segmented epidermis (also confirmed by optical inspection).This was favorable, as we expected the epidermis to be segmented as one connected component.Hence, we chose to filter with a 1D maximum filter of length 5 as a post-processing step to improve our segmentation.
To understand the quality of our segmentations, we performed an inter-rater test where two experts independently labeled 14 images of skin layers.We calculated volumetric scores of raters against each other and calculated an agreement in Dice of 0.7970 ± 0.1132 and 0.6771 ± 0.1441 in IoU.We observed that the quantitative performance of our model was on par with expert segmentations and had lower standard deviations as well.

D. Vessel Segmentation
Our neural network architecture for vessel segmentation is based on the VesNet (DeepVesselNet) architecture [32], which is a 3D fully convolutional neural network (FCN) with four fully convolutional layers.We modified the architecture by replacing 2D-crosshair convolutions with full 3D convolutions to increase performance.In addition, we added group normalization layers and increased the feature space size.We trained DeepRAP using one RTX8000 GPU.We use the following torchio augmentations: RandomFlip, RandomBlur, Random-Noise and RandomBiasField.As the final loss function (L f inal ), we used a weighted combination of the topology aware clDice [34] (L cl Dice ) and the binary cross entropy (BCE) loss function (L BC E ) to preserve vessel connectivity, where α is the weighting parameter, as in Equation 1: The networks were trained on two channels of input data (detailed in RSOM imaging and image reconstruction).We found that training on the two channels improved the segmentation performance by roughly 2% while also reducing the standard deviation (Table II).Labeling a large number of curvilinear structures such as blood vessels in 3D RSOM volumes is highly time-consuming.Therefore, we chose to follow a Transfer Learning approach for the vessel segmentation task [33].Here, we trained our neural network first on a large set of generated synthetic data samples; second, on a small set of real RSOM samples with annotated 3D vessels; and third, on so-called background samples to reduce artifacts on a set of RSOM samples without vessels.
Initially, we trained our network on synthetic arterial tree images which were generated using the method of Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

TABLE III VESSEL SEGMENTATION COMPARISON AMONG OUR MODEL WITH OTHER METHODS
Schneider et al. [44].Style transforms were applied to minimize its domain distribution shift to the original RSOM images, resulting in 30 samples with two channels and 304 × 325 × 600 pixels.To show that our synthetic dataset was indeed a valid starting set for our data, we trained networks on only synthetic data while testing real data (Table III).We observed that training purely on synthetic data leads to an acceptable segmentation performance, which can be improved with few real labeled samples (82.24 on purely synthetic data and 88.42 when refined on real data, see Table III).For the fine-tuned vessel segmentation model, 30 synthetic volumes and 15 RSOM volumes from volunteers are used for further training, 4 RSOM volumes are used for validation and 3 RSOM volumes are used for testing.The epidermis was already cut in the Z-direction, resulting in two channel data with an X, Y resolution of 333 × 171 pixels and a varying depth resolution around 400 pixels.
A frequently observed artifact in our results were layerlike reflections in the lower part of the RSOM images and above the skin surface.In response to this, four very noisy RSOM samples were selected, labeled only as "background", and used for training to improve the network's capability to distinguish between reflections and vessels.The parts of the background samples containing vascular structures were excluded beforehand so as not include false negative samples in our dataset.We extensively compared our network architecture and loss functions to other stateof-the-art methods and found that the VesNet architecture, with four convolutional layers and a combined loss function of clDice and BCE, outperformed other cost functions such as the standard Dice loss (Table III).Our method also outperforms all traditional vessel segmentation methods (e.g.Frangi, Table III).Compared to other state of the art deep learning methods, our model outperforms a simple U-Net [45], performs similar to a SwinUnet transformer [46] but is slightly outperformed by a nnUnet [31].However, we chose to use the lightweight VesNet architecture because it requires about 520 times less (30582058 vs. 58 816) parameters compared to nnUnet while still providing us with "humanlevel" performance (see Results).Using a lightweight neural network (NN) is important in a clinical setting to speed up the analysis process of the entire pipeline.Regardless, changing the segmentation NN is easily employable due to the modular setup of our pipeline.

E. Vessel Features
Based on the binary segmentations, we computed the total blood volume in mm 3 .The blood volume is defined by the sum of all segmented vessel voxels.Furthermore, we extracted the surface area to volume ratio (sa/vol), where low values describe compact shapes, and high values describe objects with large surface areas.The surface is the defined scanning field of view (4 × 2 mm 2 ).We used Pyradiomic implementations to extract these features [47].Complex vascular features were extracted from the metric graph representation, see below for details.Features were calculated by iterating over the metric graph edges or nodes.We calculated the total microvasculature length in millimeters by extracting the collective length of all edges in a volume.Using the edge weight, we distinguished small and large vascular structures, where the vessel radius r <= 2.5 pixel (30 µm) denotes a small vessel.Using this criterion, we extracted the feature of small vessel length, which is an indicator for microvasculature.The location where a vessel splits into two or more branches is a bifurcation point, a key characteristic of vascular networks.Bifurcation points were extracted from the metric graph as the number of nodes with a degree higher than two.Moreover, we extracted the average vessel radius per sample.We explored more features such as the number of loops, degree assortativity coefficient [48], and average path length, but omitted them due to low significance for subsequent use cases.

F. Metric Graph Extraction
Extracting specific features from a volumetric segmentation mask has limitations regarding a compact description of the connectivity information.Therefore, we constructed a metric graph representation (G) to achieve a compact anatomical representation preserving topological properties [49].A metric Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.graph is a one-dimensional stratified space composed of linked nodes and edges.To extract G, we followed the approach by Aanjaneya et al. [50], which is based on distinct groups of points.A pre-existing open-source project served as our codebase [50].In this work, we extended the existing 2D implementation to 3D and termed it metric graph reconstruction.As an initial step, we convert the complete segmentation mask into a point cloud.To reduce computational complexity, we constructed a skeleton based on this point cloud and passed it to our algorithm.To preserve minimal volumetric shape properties, we further extended G to hold edge weights representing the average vessel radius.Here, we determined for each element of the skeleton the closest Euclidean distance to the background and assigned this distance value to the element.The resulting modified skeleton holds the minimum vessel radius in each point.When merging groups of edge points to identify an actual edge during metric graph reconstruction, we added the mean value of all points as a weight corresponding to the edge's average vessel radius.Fig. 2 depicts the result of our modified metric graph reconstruction.We observed that metric graph reconstruction decreased the complexity of the volumetric segmentation mask by multiple orders of magnitude, while preserving major topological properties.
Due to the uniform representation of graphs, we could now easily iterate over nodes and edges and compute features quickly.

G. Post-Processing of the Graph Representation
The metric graph representation G is based on the original vessel segmentation mask (not post-processed); therefore, we discarded small unpaired structures in G with a Euclidean length < 50.Furthermore, G might contain small ending branches that do not represent the actual vascular topology caused by the skeletonization of unsmooth or spikey vessel surfaces.Thus, we "smoothed" G by removing ending branches (edges with exactly one node of degree one) with a Euclidean length smaller than 20 pixels.

H. Feature-Regions of Interest
To minimize the effects of different image sizes, varying regions of interest, and segmentation inconsistencies on the feature extraction, we applied two pre-processing steps.The most frequent misclassifications occurred slightly below the epidermis and in the lower dermal part, usually caused by socalled reflections.Thus, we restricted our feature extraction to an ROI, focusing only on the central dermal part.We estimated that the information loss due to the ROI is less severe than discarding samples that contain reflections.The ROI had the same dimensions for all volumes, allowing us to report absolute values for the features.Secondly, our feature extraction could have been affected by small false positives.Our study focused on connected vascular components and not on small, unconnected objects.Therefore, objects with a total volume smaller than 1000 pixels were removed for all features not dependent on metric graph (G).

I. Statistics
All metrics represent the mean value with standard deviations (e.g. as error bar).To assess the statistical significance between different age groups, we performed parametric tests (unpaired t-test) for normally distributed data; otherwise, nonparametric tests (Mann Whitney U test) were applied.Statistical significance was defined at P < 0.05.We applied third-order polynomial curve fitting to the data distributions of RSOM features in different age groups.

A. Quantification of Skin Layer Thickness and Vasculature Features Using Deep Learning
The superficial structure of human skin comprises of the epidermis and the dermis layers.Epidermal thickness is an important biomarker to assess the dermatological health and severity of pathologies such as psoriasis [3], [8].Therefore, we developed a CNN segmentation model inspired by our previous work [43] to segment the pathological skin layers (Fig. 1).Our U-Net achieves a high segmentation performance, with a Dice score of 84.26 ± 8.22 and IoU of 73.63 ± 11.71.These scores can be considered as "human level" performance.
We validated this via a dedicated inter-rater experiment where two experienced raters were tasked with manual segmentation of 15 identical images.The two experts achieved a Dice score of 79.70 ± 11.32 and 67.71 ± 15.62 IoU (standard deviation was calculated between the full 3D volumes), which is similar to the performance of our model (see Table I).Detailed documentation, hyper-parameter search and results can be found in the Methods section.Based on the layer segmentation shown in Fig. 2a, the average epidermal thickness in 3D was automatically computed by averaging 2D slices in the out-ofplane direction.
Furthermore, the structural and functional changes of the cutaneous microvasculature of the dermis are closely associated with changes in disease activity [3], [4], [5], [6].We developed a dedicated segmentation model for the dermal microvasculature (Fig. 1a), in which we thoroughly benchmarked and tested different architectures, layer depths, loss functions and hyperparameter configurations to achieve optimal segmentation as shown in Fig. 1b and Supplementary movie 1.Our final model is based on a DeepVesselNet architecture, which is a 3D CNN with 4-layers (Table III).We trained our model on a mixed set of images from healthy controls and patients with diabetes, using the clDice loss function [34] to preserve vessel connectivity.Our model achieved a Dice score of 88.42 ± 1.62 and 86.37 ± 6.14 using clDice.Similar to our epidermis segmentation, we again implemented an inter-rater experiment where two expert evaluators rated 16 volumes.Compared to the labels, they achieved a Dice score of 85, a performance similar to our model.Based on this segmented vessel network, various biomarkers were computed, as shown in Fig. 1.

B. Segmentation Comparison
The segmentation comparisons of various RSOM data between DeepRAP and a thresholding-based method are shown in Fig. 2. The first column (Fig. 2a) depicts the original RSOM volumes with strong reflection artifacts (marked by white arrows), which appear in the regions above the skin surface or in the lower dermis mixed with the vasculature.Since the contrast of reflection artifacts is similar to the RSOM features, the thresholding method (Fig. 2b) cannot separate the artifacts (marked by white arrows) while they are completely removed by the DeepRAP method (Supplementary movie 2).In addition, the fine microvasculature has lower contrast compared to the large RSOM features, which can be easily degraded using the thresholding segmentation (Fig. 2b).However, those fine microvasculature features are well segmented in the DeepRAP image (Fig. 2c), indicating that DeepRAP can achieve much higher segmentation accuracy with RSOM datasets containing various contrasts and artifacts compared to the thresholding method.

C. Severity Assessment of Psoriasis Skin
The epidermis thickness has been reported to be an important biomarker to diagnose and monitor the severity of psoriasis [3], [8].To investigate whether DeepRAP can differentiate between various stages of disease activity, we assessed epidermal thickness in psoriasis patients with different PASI scores ranging from 1-7 (7 being the most severe).We then applied DeepRAP to automatically segment and compute the epidermal thickness in volumetric RSOM images (Supplementary movie 3).Fig. 3a-e depict images from patients with different PASI scores with segmented epidermis masks.Using DeepRAP, we found that the PASI values are closely correlated with the epidermal thickness, i.e., the more severe the disease, the thicker the calculated epidermal thickness (Fig. 3f).Our findings are well correlated to our previous work analyzed manually by experienced RSOM operators based on MIP 2D images [8].

D. Quantification of Skin Feature Changes During Hyperemia
To demonstrate the capability of RSOM for assessing endothelial function of the cutaneous microvasculature over time, we applied DeepRAP to quantify skin feature changes during pressure-induced hyperemia.Here, we recorded 3D images of the skin structures every minute during a nineminute pressure-induced hyperemia process, consisting of a two-minute baseline measurement without cuff pressure, four minutes of cuff inflation and three minutes of cuff deflation.Cross-sectional MIP images in Figures 4a-c show that during the hyperemia process, there was a decrease in image intensity in the dermis compared to the baseline image.Following cuff deflation, we observed that vascular structures fully recover, with dilated vessel diameters and more cutaneous vessels (white arrows in Fig. 4c).It is worth noting that the skin microvasculature was accurately segmented despite strong variations in image contrast during the hyperemia process.In addition, the microvascular changes during hyperemia were characterized by computing the mean RSOM image intensity (Fig. 4d).Next, we applied DeepRAP to measure epidermal thickness and various vascular features including the total blood volume (Fig. 4e), total vessel number (Fig. 4f), total vessel lengths (Fig. 4g) and small (diameter less than 60 µm, Fig. 4h) and large vessel (diameter more than 60 µm, Fig. 4i) lengths.The response pattern of these vascular biomarkers was Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.The normalized image intensity profiles of the RSOM images acquired during the hyperemia process.Skin features computed from the segmented images include: (e) total blood volume, (f) total vessel number, (g) total vessel length, (h) total length of small vessels (with diameter less than 60 µm), (i) total length of large vessels (with diameter more than 60 µm), (j) diameter change of the labelled vessel (marked by the white dashed line in a), (k) epidermal thickness; (l) profiles of the blood flow, oxygen saturation (SO2), and partial blood volume (rHb).Scale bar: 500 µm.
similar to the image intensity profile (Fig. 4d).In addition, the changing diameter of a specific vessel (marked by the white dashed line in Fig. 4a) during the hyperemia process was well visualized and quantified in Fig. 4j.Arterial occlusion did not affect epidermal thickness as it remained relatively constant (Fig. 4k).To compare our results, we used a commercial Laser Doppler flowmetry and tissue spectrometry setup to measure the blood flow (Flow), oxygen saturation (SO 2 ) and partial blood volume (rHb) during post-occlusive reactive hyperemia (PORH).The changes in flow and SO 2 were similar in trend to the skin biomarkers computed by DeepRAP (Fig. 4l).These change profiles of the individual parameters can be used to quantify skin vessel function, for example, the endothelial function of the small and large vessels respectively.

E. Assessment of Aging on Skin Features
Previous studies have reported that skin morphological and microvasculature features can be affected by aging [13], [51], [52], [53], [54].To further validate the utility of our deep learning tool, we studied skin feature changes resulting from aging, applying DeepRAP to analyze RSOM datasets acquired from healthy volunteers with different ages.Fig. 5a-c shows cross-sectional MIP RSOM images and corresponding dermal vessels from three groups: 1. young age volunteers (I, n=24, 29.5 ± 3.5 years, in the range of 20 to 35 years), 2. middle age volunteers (II, n= 28, 41.1 ± 4.1 years, in the range of 36 to 50 years) and 3. old age volunteers (III, n=23, 59.5 ± 4.9 years in the range of 51 to 65 years).Using DeepRAP, we found that the density of the dermal vascular network was highest in young volunteers (Fig. 5a) and decreased with age, with lower densities found in volunteers from the middle age group (Fig. 5b) and the lowest densities found in the old age volunteers (Fig. 5c).Vessel lengths, diameter and the number of bifurcation points are commonly computed to quantify the vascular architecture.Hence, we used DeepRAP to segment the RSOM volume and compute these vascular features and the epidermal thickness.Quantitative comparisons of the magnitude of different RSOM features are presented in Fig. 5d-5h, while the distributions of these features with the increment of age is shown in Fig. 5i.As for the total vessel length (Fig. 5d, length normalized to the size of the scanning field of view), we found a mean value of 21.43 ± 4.48 mm in the young age (I) group versus 16.52 ± 5.62 mm in the middle age (II) group (p = 0.0022) with 23% difference of the mean value, while the mean total vessel length significantly decreased to 10.19 ± 2.39 mm in the old age (III) group, a 38% reduction compared to the middle age (II) group (p < 0.001).To investigate the effects of age on different vessel sizes, we further separated the total vessel length into the small (diameter < 60 µm) and large (diameter ≥ 60 µm) vessel lengths.The small and large vessel lengths both decreased significantly from the young age group to the middle age groups (for small vessels, 5.25 ± 1.84 mm (I) vs. 3.58 ± 1.49 mm (II), with p = 0.0012; and for large vessels, 17.17 ± 3.83 mm (I) vs. 12.84 ± 4.89 mm (II) with p = 0.0288).We found that the average small vessel length (Fig. 5e) was approximately 56% lower in the old age group (III) than in the middle age (II) group (1.57± 0.46 mm vs. 3.58 ± 1.49 mm with p < 0.0001), whereas the mean large vessel length (Fig. 5f) was about 49% lower in the old age group (III) compared to the middle age (II) group (8.75 ± 2.34 mm versus 17.17 ± 3.83 mm, p = 0.0019).This finding suggests that the systemic impacts of aging on the dermal vasculature are more prominent in small vessels than in larger vessels.The total blood volume in the DR layer (i.e., the sum of the total vessel voxels) was markedly different between the young and middle age groups (Fig. 5g, 0.043 ± 0.0074 mm 3 versus 0.034 ± 0.0093 µm 3 , p < 0.0008, 21% mean value difference), while the difference is also larger between the old age and the middle age groups (0.02 ± 0.0067.mm 3 versus 0.034 ± 0.0093 mm 3 , p < 0.0001, 40% mean value difference).The impairment of aging on the dermal microvasculature is more significant between the middle age to old age groups compared to the young age to middle age groups.Analysis of the epidermal thickness (Fig. 5h) showed no obvious difference among the young age (109.80 ± 12.37 µm), middle age (108.76 ± 9.94 µm) and old age groups (109.69 ± 14.35 µm).

IV. DISCUSSION
DeepRAP is a pipeline for automatic segmentation and biomarker extraction of the skin layers and dermal vasculature, facilitating a new avenue for disease characterization.We show excellent, human-level performance with computed results which correlate well with clinical psoriasis severity scores.We also showed that changes in skin morphology and vasculature during PORH have similar response patterns of blood flow and oxygen saturation as measured by commercial setups.These experiments demonstrate the efficacy of DeepRAP for comprehensive assessment of microvascular endothelial function facilitating clinical application [55].In addition, five anatomical and vascular features are extracted by DeepRAP and studied in relation to aging.We found that the small vessels in the upper dermis were most prominently impaired with aging.
In previous RSOM studies, conventional analysis often relied on manual or thresholding-based segmentation methods to process 2D MIP images, which could induce significant variations, annotator bias or loss of the 3D geometrical information.This method is error-prone and labor-intensive when used to analyze large clinical RSOM Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
studies.DeepRAP is the first deep learning-based method designed to automatically segment RSOM images in 3D and compute various skin features.We also demonstrated that DeepRAP achieves much better segmentation performance compared to conventional thresholding-based methods (Fig. 2) and is more robust towards varying SNR and artifacts.In a unique validation experiment, we applied DeepRAP to a timeseries sequence of RSOM volumes recorded over occlusion-induced hyperemia from 10 healthy volunteers.We observe how the biomarkers decrease and recover during the occlusion and release process, demonstrating accurate performance and reproducibility in a highly challenging validation task which exhibits marked contrast variations.Previous studies report an urgent need to develop a safe and ideally noninvasive method of microvascular function assessment [55] which we address in this study.Hence, assessment of the skin micro-vasculature endothelial function is important for early disease detection and therapy monitoring, e.g., for diabetes and cardiovascular diseases [56], [57].Our DeepRAP pipeline can accurately segment and quantify skin features from a series of RSOM images, despite significant changes in image contrast during hyperemia test.This allows RSOM to assess microvasculature function by accurately quantifying vessel changes in response to stimuli, which could further promote the clinical applications of RSOM.Furthermore, DeepRAP was able to analyze large RSOM datasets from healthy and disease conditions, with the capacity to automatically measure skin morphology features in 3D for disease monitoring in the relevant anatomical regions of interest.
Using DeepRAP, we showed that vessel lengths (of varying hierarchy), vessel number and total blood volume, were reduced in the middle age group compared to the young age group and significantly decreased in the old age group compared to the middle age group which is in line with previous findings [13], [51], [52], [53], [54].For example, it was reported that the superficial skin microvasculature assessed by video capillaroscopy or histology was significantly reduced in older volunteers compared to young volunteers.Moreover, aging impaired the small vessels more significantly than the large vessels, which has not been reported in previous studies.
DeepRAP is based on a transfer learning approach, where we pre-trained the CNN on synthetic data [44] and refined it on a small labelled dataset of 17.75% of the synthetic dataset.Thus, our method might generalize well to different types of imaging data (such as other optoacoustic imaging systems or other tissue structures), as only a small, labelled dataset is needed to adjust our pre-trained network.It has been shown that for the segmentation of tubular structures such as vessels, the optimization of voxel overlap alone (e.g., Dice) is not sufficient [34], [58].Hence, we aim to improve the topological faithfulness of our method by employing a topology-aware loss function, namely clDice.However, this is only one of many approaches for improving topology.Frequently used methods are based on post-processing [59], tree shape priors for vessel segmentation [60], or persistent homology [37].Approaches based on persistent homology come with the strictest theoretical guarantees.Nevertheless, in our use case, we chose clDice because it is more robust towards image noise, as present in RSOM data, which presents challenges when modeling topology in structures used in persistent homology approaches, such as Betti numbers.Furthermore, existing approaches address vessel analysis directly at a graph level (segmentation, in our case), allowing the computation of total blood volume (Fig. 5), an important biomarker.From a limitation perspective, it is important to be aware that all supervised learning methods are heavily dependent on data and label quality.Secondly, while transfer learning helps with the generalizability of methods, this generality has limits, and the performance of the model should always be evaluated on a target domain test set.Finally, training topological methods requires more computational resources than training basic networks with, for example, Dice loss.
In conclusion, DeepRAP is a scalable, modular and automated machine learning-based method that can be used to analyze and quantify skin morphological and functional features from 3D RSOM datasets.We envisage that our method will be employed to promote the applications of RSOM imaging for quantitative clinical studies and diagnostics.

Fig. 1 .
Fig. 1.DeepRAP processing diagram.(a) the reconstructed high frequency (HF) and low frequency (LF) images are preprocessed separately as the input data; the raster-scan optoacoustic mesoscopy (RSOM) image volumes are segmented using a U-Net to precisely delineate the epidermis and dermis layers and the dermal vasculature is segmented by a VesNet; vascular graph is extracted, and various skin biomarkers are computed.(b) Segmentation results of one RSOM volume; 3D renderings of the original RSOM volume (left); results of our deep learning layer and vessel segmentation (middle), where the segmented epidermis (EP) is marked in green, and the segmented dermal (DR) vessels are in red; dermal vessels are color coded according to their diameter (right).The bottom row shows the maximum intensity projection images from the top view corresponding to the dermal vessels.HF, high frequency; LF, low frequency.Scale bar: 500 µm.

Fig. 2 .
Fig.2.Multiple segmentation comparisons on raster-scan optoacoustic mesoscopy (RSOM) images between our DeepRAP and a thresholdingbased approach.Using DeepRAP, the segmented epidermis (EP) is marked in green, and the segmented dermal (DR) vessels are in red.A threshold-based approach cannot distinguish these.(a) DeepRAP accurately segments the epidermis and dermis microvascular layers of RSOM volumes with artifacts (white arrows) mixed with the epidermis layer while these artifacts are left in the thresholding image; (b) DeepRAP segments the microvasculature in the high frequency content (yellow arrows) while these vascular features are not correctly segmented using the thresholding approach; (c) DeepRAP accurately segments the dermal microvasculature by removing artifacts mixed with dermal microvasculature while they are segmented as vessels in the thresholding image (white arrows).

Fig. 3 .
Fig. 3. Epidermis segmentation of psoriasis skin.(a)-(e) Cross-sectional raster-scan optoacoustic mesoscopy (RSOM) images of psoriasis skin lesions with various Psoriasis Area Severity Index (PASI) scores, where the epidermal thickness is segmented and marked by the grey areas.(f) Correlation between the PASI score and the epidermal thickness.Scale bar: 500 µm.

Fig. 4 .
Fig. 4. Quantification of skin features during hyperemia.3D RSOM images were recorded at the forearm of a healthy volunteer every minute during a nine-minute arterial occlusion process.The original and segmented RSOM images by DeepRAP at (a) the baseline, (b) cuff inflation, and (c) deflation points are shown.White arrows indicate new vessels appearing during the hyperemia process.The epidermis is segmented and marked by the grey areas, while the dermal vessels are marked in blue color and overlaid with the original vessel structures.(d)The normalized image intensity profiles of the RSOM images acquired during the hyperemia process.Skin features computed from the segmented images include: (e) total blood volume, (f) total vessel number, (g) total vessel length, (h) total length of small vessels (with diameter less than 60 µm), (i) total length of large vessels (with diameter more than 60 µm), (j) diameter change of the labelled vessel (marked by the white dashed line in a), (k) epidermal thickness; (l) profiles of the blood flow, oxygen saturation (SO2), and partial blood volume (rHb).Scale bar: 500 µm.

Fig. 5 .
Fig. 5. Comparisons of skin features among healthy volunteers in three age groups.The healthy volunteers were divided into three groups based on their ages: I (n=24, 29.5 ± 3.5 years), II (n= 28, 41.1 ± 4.1 years) and III (n=23, 59.5 ± 4.9 years).Three cross-sectional raster-scan optoacoustic mesoscopy (RSOM) images of each group and corresponding dermal vessels from the top-view are shown in (a-c), where the epidermis (EP) and dermal (DR) vessels were segmented by the DeepRAP.The epidermis is segmented and marked by the grey areas, while the dermal vessels are marked in blue color and overlaid with the original vessel structures.(d-h).Skin features were computed and compared among the three groups including: (d) total vessel length, (e) small vessel length (vessels with diameter < 40 µm), (f) large vessel length (vessels with diameter ≥ 40 µm), (g) total blood volume and (h) epidermal thickness.(i) the distributions of these five skin features with increment of volunteer aging.ns: not significant.Scale bar 500 µm.