Estimating quantitative physiological and morphological tissue parameters of murine tumor models using hyperspectral imaging and optical profilometry

Understanding tumors and their microenvironment are essential for successful and accurate disease diagnosis. Tissue physiology and morphology are altered in tumors compared to healthy tissues, and there is a need to monitor tumors and their surrounding tissues, including blood vessels, non‐invasively. This preliminary study utilizes a multimodal optical imaging system combining hyperspectral imaging (HSI) and three‐dimensional (3D) optical profilometry (OP) to capture hyperspectral images and surface shapes of subcutaneously grown murine tumor models. Hyperspectral images are corrected with 3D OP data and analyzed using the inverse‐adding doubling (IAD) method to extract tissue properties such as melanin volume fraction and oxygenation. Blood vessels are segmented using the B‐COSFIRE algorithm from oxygenation maps. From 3D OP data, tumor volumes are calculated and compared to manual measurements using a vernier caliper. Results show that tumors can be distinguished from healthy tissue based on most extracted tissue parameters ( p<0.05 ). Furthermore, blood oxygenation is 50% higher within the blood vessels than in the surrounding tissue, and tumor volumes calculated using 3D OP agree within 26% with manual measurements using a vernier caliper. Results suggest that combining HSI and OP could provide relevant quantitative information about tumors and improve the disease diagnosis.

tumors compared to healthy tissues, and there is a need to monitor tumors and their surrounding tissues, including blood vessels, non-invasively. This preliminary study utilizes a multimodal optical imaging system combining hyperspectral imaging (HSI) and threedimensional (3D) optical profilometry (OP) to capture hyperspectral images and surface shapes of subcutaneously grown murine tumor models. Hyperspectral images are corrected with 3D OP data and analyzed using the inverseadding doubling (IAD) method to extract tissue properties such as melanin volume fraction and oxygenation. Blood vessels are segmented using the B-COSFIRE algorithm from oxygenation maps. From 3D OP data, tumor volumes are calculated and compared to manual measurements using a vernier caliper. Results show that tumors can be distinguished from healthy tissue based on most extracted tissue parameters (p < 0:05). Furthermore, blood oxygenation is 50% higher within the blood vessels than in the surrounding tissue, and tumor volumes calculated using 3D OP agree within 26% with manual measurements using a vernier caliper. Results suggest that combining HSI and OP could provide relevant quantitative information about tumors and improve the disease diagnosis.

| INTRODUCTION
Understanding tumors and tumor microenvironment are essential for successful diagnostic and treatment approaches. Hanahan and Weinberg [1] have initially postulated the six hallmarks of cancer, that is, distinctive biological capabilities that enable tumor growth and metastatic dissemination. The fundamental trait of cancer cells is their ability to sustain chronic proliferation (rapid reproduction), evade growth suppressors, resist cell death and activate invasion and metastasis. Recently, the six hallmarks were expanded, and among the new hallmarks of cancer are inflammation and metabolic specificities [2]. Cancer cells require large amounts of nutrients and oxygen and therefore induce angiogenesis, the formation of new blood vessels, to meet the increased demand. Unlike in normal tissue, the blood vessels in tumor tissue are immature, tortuous and highly disorganized, resulting in hyperpermeability (leakage), low perfusion, hypoxia (deprivation of oxygen) and decreased blood flow within the tumor [3].
Tumor vasculature presents a complex hurdle for achieving therapeutic success across various cancer treatment approaches. The aberrant vasculature presents a physical barrier to drug treatment and promotes aggressive and invasive tumor phenotypes because of the tumor's hypoxic environment. The abnormal blood circulation leads to insufficient perfusion and drug delivery and creates pockets of hypoxia, rendering solid tumors more resistant to radiotherapy (RT), chemotherapy and immunotherapy [3][4][5]. The undesired therapy side effects such as decreased radiosensitivity due to hypoxia in radiotherapy highlight the pressing need for studies of tumors and their vasculature.
Optical imaging techniques are sensitive to light absorption and scattering changes and can discriminate healthy tissue from tumor tissue based on the intrinsic tissue contrasts [6,7]. Hyperspectral imaging (HSI) is a promising non-invasive and contactless optical imaging technique that captures spatial and spectral information of the imaged tissue in a hyperspectral image (hypercube), commonly in the visible (VIS) and near-infrared (NIR) spectral bands [8]. The imaging system can operate in different modes to capture the reflection or transmission spectra of the observed biological tissue. The resulting hypercubes contain spectral information in each pixel of an image. Following the hyperspectral image analysis, diagnostic information about tissue physiology (e.g., melanin volume fraction or tissue oxygenation), pathology (e.g., increased metabolic activity) and morphology (e.g., tissue components size) is obtained. [8] Based on the extracted parameters, HSI can be used to differentiate between healthy and tumor tissue in vivo or extract blood vessel maps [8][9][10].
In recent years, HSI has been widely used to image human and animal tumors in brain [8,9,[11][12][13][14][15][16][17][18][19][20][21][22], breast [8,9,11,[23][24][25][26][27][28], colon [8,9,11,[29][30][31][32][33], prostate [8,9,11,34,35] and skin [8,9,11,[36][37][38][39][40][41][42]. The main shortcoming of most studies utilizing HSI is that they analyze images of histopathological slides to detect tumor cells ex vivo. Although this approach fits well in the current standard of care where biopsy is taken, and tissue histology is inspected, this approach has many limitations. First, such a procedure is invasive as tissue samples must be taken. Moreover, the biopsied tissue's optical properties differ from the intact in vivo tissue, changing over time due to dehydration, the absence of blood and tissue fixation [43,44]. Thus, many tissue parameters cannot be extracted or are altered compared to the intact tissue. Finally, the progression of angiogenesis and blood vessel growth within the tumor and its surroundings cannot be monitored due to the invasive nature of the biopsy. On the other hand, macroscopic in vivo hyperspectral imaging could image a large area (a few cm 2 ) of tumor tissue in a non-invasive and contactless fashion, leaving the tissue completely intact during and after the procedure. Additionally, the time to process the information and establish a diagnosis is significantly reduced by eliminating the need to prepare histopathological samples.
HSI can be complemented with other non-invasive optical techniques to provide additional information about the biological tissue. A three-dimensional optical profilometry (OP) provides the tissue surface shape [45] that can be used to accurately calculate tumor volumes and eliminate the need to measure tumor sizes with a vernier caliper, a standard procedure [46,47] in assessing tumor volumes. Although optical profilometry is a widespread technique in industry and many scientific fields, it is still not common in biomedical optics to extract information about tumors or other biological tissue. Using the surface shape information, the tumor topology, such as size and shape, can be assessed, and tumor texture can be extracted based on the roughness of its surface. Additional information could help differentiate between tumors and healthy tissue or different tumor types. For example, Handels et al. [48] showed that laser profilometry could improve non-invasive melanoma diagnosis in human patients. Norhaimi et al. [49] demonstrated that a fringe projection imaging modality could detect changes in breast surface caused by the size variation of the tumor in affected women. In another study, Pavlovčič et al. [50] applied a compact and rotational laser profilometer to measure the volume of skin wounds accurately.
The main aim of this preliminary study is to noninvasively determine quantitative physiological and morphological tissue parameters of four different murine tumor models in vivo and compare them to healthy tissue. Doing so could help discriminate the tumors from the healthy tissue and differentiate various tumor models. Therefore, this study utilized macroscopic hyperspectral imaging to record spatially resolved spectral images of tumors and applied a GPU-accelerated twolayer inverse adding-doubling (IAD) method to determine tissue properties from hyperspectral images. What is more, tumor growth is closely connected to the process of angiogenesis and blood vessel growth within the tumor and its surroundings to ensure a sufficient supply of oxygen and nutrients. Thus, one of the main goals was also to automatically segment the blood vessel maps using the B-COSFIRE [51][52][53] algorithm to determine blood oxygenation within the blood vessels. Finally, a standard operating procedure in studying murine tumor models involves manually measuring the tumor dimensions using a vernier caliper to determine tumor volumes. Hence, another objective of this study was to obtain the surface shape of the imaged tissue to determine tumor volumes using the optical profilometry module integrated with the HSI system and compare them with manual measurements. This part aims to facilitate and increase the accuracy of the tumor size determination in animal tumor models.

| Imaging system
This study used a custom-built multimodal optical imaging system combining a hyperspectral imaging (HSI) module and a three-dimensional optical profilometry (OP) module presented in Figure 1. HSI was integrated with an optical profilometry module to obtain the surface shape of the imaged tissue and apply the curvature and height corrections [45] of hyperspectral images to compensate for the signal loss in areas with high inclination angles and distance variation. HSI is a push-broom (linescanning) hyperspectral imaging device consisting of a monochrome CMOS camera (Blackfly S, BFS-U3-51S5M-C, FLIR, Canada), an imaging spectrograph (ImSpector V10E, Specim, Finland), two motorized translation stages (8MT195, Standa, Lithuania), custom LED illumination panel spanning the spectral range of 400-1000 nm and crossed polarizers to minimize specular reflection. A 50 mm objective (Xenoplan 2.8/50-0902, Schneider-Kreuznach, Germany) is used to image a small field of view (FOV) suitable for murine tumor models. The 3D optical profilometry module was based on the triangulation method and comprised a laser projector (FLEXPOINT, 30 mW, 405 nm, LASER COMPONENTS, Germany) and a monochrome camera (Flea3, FL3-U3-13Y3M-C, FLIR, Canada) with a 16 mm lens and a bandpass filter.
The HSI and OP modules were aligned using a reference object of known geometry resulting in the multimodal image misalignment <0.1 mm. The spatial and spectral resolution of the HSI module was 100 μm and 2.9 nm, respectively. The accuracy of a 3D surface captured by the OP module in the X, Y and Z directions was 100, 100 and 50 μm, respectively.

| Experimental procedure
Seven BALB/c (BALB/cAnNCrl, Charles Rivers) and one C57Bl/6 (C57Bl/6NCrl, Charles Rivers) 8-10 weeks old female mice were included in the study to image subcutaneously grown tumors. Specifically, two mice with a 4T1 murine mammary carcinoma (American Type Culture Collection, ATCC), one mouse with a B16-F10 murine melanoma (ATCC), three with a CT26 murine colon carcinoma (ATCC) and two with a TS/A [54] mouse mammary adenocarcinoma (authenticated by CellCheck, mouse short tandem repeat [STR] profile and interspecies contamination test) were imaged in the preclinical setting at the Department of Experimental Oncology, Institute of Oncology Ljubljana. Mice were kept in a specific pathogen-free environment with a 12-h light-dark cycle at 20 C-24 C with F I G U R E 1 A schematic of the multimodal optical imaging system combining HSI and OP 55% ± 10% relative humidity and food and water provided ad libitum. The experiments were approved by the Ministry of Agriculture, Forestry and Food of the Republic of Slovenia (permission no. U34401-3/2022/11).
Before imaging, mice were anesthetized with an anesthetic solution (ketamine, 125 mg/kg; xylazine 12.5 mg/kg; acepromazine 2.5 mg/kg) injected intraperitoneally with volumes adjusted to the weight of the animal. The back of the mouse was then shaved and depilated with a depilatory cream to expose the bare skin and avoid excess light scattering on the white hair.

| Image preprocessing
Raw hyperspectral images were normalized to calculate reflectance using [8]: where I ref is the sample reflectance, I raw is the raw hyperspectral pixel intensity, I dark is the dark current and I white is the white standard intensity (Spectralon, Labsphere Inc., New Hampton). Normalized hyperspectral images were then corrected for the signal loss due to the changing inclination of the object surface and distortions resulting from the varying distance between the object and the camera. Specifically, 3D OP data were analyzed to obtain the surface shape profile, and height and surface curvature corrections were applied as described by Rogelj et al. [45] to flatten the hyperspectral intensity profiles at different surface inclinations. Corrected hyperspectral images were spectrally cropped and binned to 430-750 nm spectral range and 5 nm spectral resolution and spatially binned tenfold in both directions to facilitate the image analysis. Ultimately, the spectral angle mapper [8,55] (SAM) was used for hyperspectral image segmentation to remove the background and obtain the segmentation masks of murine tumor models. SAM calculates the spectral similarity of two spectra by treating them as vectors in space with dimensionality equal to the number of wavelengths and calculating the angle between them [55]: where s 1 ! and s 2 ! are the two spectra of interest. As a threshold cos θ ¼ 0:2 was selected, resulting in precise segmentation of tumors from the background.

| Image analysis by the inverse adding-doubling algorithm
In order to extract physiological and morphological information from the corrected reflectance images, the biological tissue (e.g., murine skin) is modeled, taking into account its geometry, absorption and scattering properties, and light propagation within the tissue must be simulated. [7] The most common techniques to extract tissue properties are the inverse Monte Carlo (IMC), the gold standard, the inverse diffusion approximation (IDA) and the inverse adding-doubling (IAD) [6]. IMC method is the most widely used, flexible and accurate, but requires high computational power and is generally time-consuming. IDA is simple and rapid but considerably inaccurate. Finally, IAD is much faster and as accurate as IMC but only applies to layered tissue [56]. Using these techniques, tissue properties are extracted iteratively until the calculated values of reflectance or transmittance do not match the measured values [6]. A GPU-accelerated inverse adding-doubling [56] (IAD) algorithm was implemented in MATLAB R2020b (Mathworks, MA), allowing rapid and accurate light propagation simulation in the murine skin model. A twolayer skin model presented in Figure 2 was used, consisting of the epidermis and semi-infinite dermis with 11 tissue properties thoroughly described in our recent publication [57]. Briefly, the absorption coefficient of the epidermis was calculated as [7]: F I G U R E 2 A two-layer murine skin model consisting of the epidermis and semi-infinite dermis with 11 tissue parameters where f m is the volume fraction of melanin, μ a,m is the melanin absorption coefficient and μ a,base is the baseline absorption of bloodless skin. The absorption coefficient of the semi-infinite dermis was defined as [7]: where f Hb and f HbO 2 are volume fractions of deoxy-and oxyhemoglobin, μ a,Hb and μ a,HbO 2 are corresponding absorption coefficients, and f brub and μ a,brub are concentration and absorption coefficient of bilirubin, respectively. Moreover, f CO and f COO 2 are concentrations of reduced and oxidized cytochrome C oxidase, whereas μ a,CO and μ a,COO 2 are associated absorption coefficients. These two parameters were included in the model to improve the fitting for wavelengths above 650 nm. The reduced scattering coefficient was calculated as [7]: Here, a is the reduced scattering coefficient at 500 nm, b is an exponential parameter related to the size of the Mie scatterers and f Ray represents the fraction of Rayleigh scattered light. In Figure 2, d e represents the epidermis thickness, and d d represents the actual dermis thickness, restricted in the IAD algorithm to values an order of magnitude thicker than the penetration depth of light in the 400-1000 nm spectral range (a few cm compared to a few mm). Since d d is so thick (the average value estimated by the IAD algorithm is 5.62 ± 0.69 cm) that almost no light passes through it can be considered semi-infinite. Levenberg-Marquardt (LM) algorithm was adopted to GPU for least-squares fitting of measured reflection spectra, and the maximum number of iterations was limited to 200. Incoming and outgoing light was divided into 20 fluxes to ensure accuracy comparable to MCML. [58] Fitting was performed on a computer comprising an Nvidia Titan Xp graphics card with 12 GB RAM, AMD Ryzen 7 1700X processor and 16 GB RAM. The batch size (simultaneous number of GPU threads) was 1000. The average number of iterations was 25 AE 14, and the average fitting time per spectrum was 0:22 AE 0:57 s.

| Blood vessel segmentation
Automatic blood vessel segmentation was performed to obtain maps and skeletons of blood vessels from blood parameter maps extracted from hyperspectral images.
First, blood oxygenation was calculated from deoxy-and oxyhemoglobin volume fractions: Then, StO 2 parameter map was thresholded and enhanced using contrast limited adaptive histogram equalization (CLAHE) to improve the contrast of blood vessels. Images were also padded to minimize the edge artifacts [10]. Finally, blood vessels were segmented using an automated B-COSFIRE [51][52][53] algorithm, which detects elongated features such as lines (i.e., blood vessels) using the identification of maximal values of convolution of the original image with a linearly arranged Difference of Gaussian (DoG) functions. Images resulting from the B-COSFIRE segmentation were thresholded to discard the parts of an image (with a response lower than 30) corresponding to noise and algorithm artifacts. Segmented binary blood vessel masks were post-processed using dilation, erosion and Gaussian smoothing to remove the artifacts and improve the segmentation. Ultimately, vessel skeleton maps were extracted from the segmented blood vessel maps.

| Tumor volume estimation
Tumor volumes were estimated using a standard procedure involving manual measurements performed at the Department of Experimental Oncology, where three perpendicular axes (A, B and C) are measured with a vernier caliper (Mitutoyo, Japan) and the tumor volume is approximated with an ellipsoid: These results were compared to the tumor volumes calculated from the three-dimensional point clouds obtained with the profilometry module of the multimodal imaging system. Two approaches were used: in the first approach, the perpendicular axes of the tumor were measured from the 3D point clouds, and tumor volumes were calculated using Equation (7); in another approach, the tumor mask from SAM segmentation was applied to 3D point cloud to locate the tumor. Then, we summed the number of voxels under the mask belonging to the tumor mass. The tumor volume was then estimated from a total number of tumor voxels and a known volume of a single voxel (V 1 ) as: 3 | RESULTS

| Image analysis and tissue parameters extraction
In this section, the results of hyperspectral image analysis of eight mice with subcutaneously grown tumors are presented. Figure 3 shows the in vivo measured reflectance skin spectra at 430-750 nm of four different tumors (in red) and surrounding healthy tissue (in blue) extracted from hyperspectral images. The following tumor models are presented: (A) 4T1, (B) B16-F10, (C) CT26 and (D) TS/A. The figure shows that healthy tissues feature either a double camel hump (oxygenated tissue) or a dip (relatively deoxygenated tissue) in the 550-600 nm spectral region. However, this does not apply to the healthy tissue surrounding the B16-F10 tumor, which is deoxygenated and has much lower reflectance values than other healthy tissues for wavelengths larger than 650 nm. As for tumor models, the reflectance spectra do not have the described characteristic shape. Generally, they have lower reflectance values, indicating higher absorption in the tumor tissue at low wavelengths, especially in the B16-F10 murine melanoma model, where high melanin volume fraction yields reflectance values close to 0 ( Figure 3B). On the other hand, the spectra of the CT26 tumor model are similar to spectra of healthy tissue, yet they are slightly more oxygenated with a higher blood volume fraction, as indicated by a more pronounced double camel hump and lower reflectance for wavelengths lower than 650 nm. It is also interesting to note that spectra of 4T1 and TS/A tumor models are very similar, indicating a similar composition of both tumor models.
In general, the spectra of the tumors could be easily discriminated from the spectra of the healthy tissue allowing the classification of tumorous and healthy tissue. The spectra also show that some tumor models, such as TS/A, are more homogeneous than the others (e.g., 4T1), deducing from the high standard deviations.
The measured reflectance skin spectra were fitted using the IAD algorithm to extract tissue parameters from the hyperspectral images. Figure 4 shows the measured and fitted spectra of a female BALB/c mouse with a subcutaneously grown 4T1 mammary carcinoma. The presented spectra are from the less oxygenated ( Figure 4A), the more oxygenated area of the healthy tissue ( Figure 4B), and the tumor tissue ( Figure 4C). The difference in oxygenation is notable in the 550-600 nm spectral region, where a double camel hump is more pronounced in Figure 4B than in Figure 4A. The tumor spectrum is significantly different from the healthy tissue spectrum. An excellent match between fitted and measured spectra is obtained (R 2 > 0:99 for the presented spectra), apart from the minor deviations for wavelengths around 400 nm and above 650 nm. The latter could be attributed to higher noise levels in these spectral regions due to the LED light source. In these spectral regions, the illumination intensity is reduced [59].  Figure 6 presents maps of melanin volume fraction f m , deoxyhemoglobin f Hb and oxyhemoglobin f HbO 2 , extracted from the hyperspectral image of a female BALB/c ( Figure 5). Notably, the melanin and deoxyhemoglobin volume fractions are much higher within the tumor tissue than in the neighboring healthy skin tissue, as are some areas of oxyhemoglobin. Blood vessel architecture is distinctly visible from the colormaps of deoxyand oxyhemoglobin (marked with black arrowheads).
The mean value of f m within the tumor and in the healthy tissue was 3.24% ± 2.51% and 0.49% ± 0.07%, respectively. The respective mean values of blood oxygenation in the tumor area and healthy tissue were 13.47% ± 18.51% and 25.52% ± 14.49%. These results suggest that the 4T1 mammary carcinoma in a murine tumor model has a significantly higher melanin volume fraction and is considerably less oxygenated than healthy skin tissue surrounding the tumor. However, high f m could be due to an interplay between melanin and necrosis, as discussed later. Although the maps are generally smooth, they appear grainy and noisy in the tumor due to the higher heterogeneity of the tumor tissue compared to the healthy tissue.
In further analysis, all tumors and healthy tissues were considered together to deduce some general tumor properties that could differentiate them from healthy tissue. Shown in Figure 7 are the box plots for average tissue parameters for all tumor models and healthy tissues. Interestingly, values of physiological properties such as f m , f Hb , f HbO 2 and StO 2 are all higher in tumors than healthy tissues, suggesting that the investigated tumors generally contain more melanin and are more oxygenated than healthy tissue. The difference in melanin volume fraction is especially relevant in the B16-F10 murine melanoma, which is six times higher within the F I G U R E 4 In vivo measured (blue line) and fitted (red line) reflectance skin spectra of a BALB/c mouse (subject 2) with a subcutaneously grown 4T1 mammary carcinoma: A, less oxygenated tissue; B, more oxygenated tissue; C, tumor tissue F I G U R E 5 Color image extracted from the hyperspectral image of a female BALB/c mouse (subject 2) with a subcutaneously grown 4T1 mammary carcinoma melanoma than in the surrounding tissue. Tissue oxygenation was higher in all examined tumor models except for a 4T1 mammary carcinoma presented in Figure 6, being an outlier for StO 2 in Figure 7D. In this tumor model, there is an area of low oxygenation in the necrotized tissue, vastly contributing to overall low oxygenation within the tumor. It could also be attributed to pockets of hypoxia within the tumor formed due to oxygen requirements exceeding the oxygen supply. We have also found that total hemoglobin content was, on average, four times higher in tumors than in healthy tissue. As for the scattering properties connected with tissue morphology, the scattering coefficient and the scattering power were higher in the healthy tissue than in tumors, suggesting a more homogenous composition of the tumors.
We performed the Mann-Whitney U-test to test if there are statistically significant differences between tumors and healthy tissue based on the extracted tissue parameters. Briefly, the test checks the null hypothesis that data in two groups are samples from continuous distributions with equal medians against the alternative hypothesis that they are not. The test assumes that the two samples are independent. For all parameters presented in Figure 7, as well as total hemoglobin content, there are statistically significant differences (p < 0:05) between tumors and healthy tissues, except for the scattering amplitude a (p ¼ 0:13). Based on these results, the spectral analysis clearly shows F I G U R E 6 Maps of A, f m ; B, f Hb and C, f HbO2 extracted from the hyperspectral image of a female BALB/c mouse (subject 2) with the subcutaneously grown 4T1 mammary carcinoma. Black arrowheads point at the visible blood vessels, and red arrowheads point at the necrotized tumor area. Color bar for f Hb is in logarithmic scale for the optimal display of blood vessels F I G U R E 7 Box plots of extracted tissue parameters for all healthy tissues and tumor models: A, f m ; B, f Hb ; C, f HbO2 ; D, StO 2 ; E, a and F, b. Tissue oxygenation, StO 2 , was calculated from Equation (6). HT, healthy tissues; T, tumors. On each blue box, the central mark represents the median value, and the bottom and top edges indicate 25th and 75th percentiles, respectively. The whiskers extend to the most extreme data points not considered outliers, and the outliers are plotted individually using the red plus (+) marker symbol that the tumors have significantly different physiological and morphological properties than healthy skin tissue. Therefore, spectral imaging enables the discrimination between them.

| Blood vessel segmentation
Herein we present the results of blood vessel segmentation from oxygenation maps extracted from hyperspectral images of murine tumor models. As described in Section 2.5, tissue oxygenation was first calculated using Equation (6), followed by image contrast enhancement using CLAHE and blood vessel segmentation using the B-COSFIRE algorithm. Shown in Figure 8 are the results of blood vessel segmentation from calculated StO 2 parameter maps. Specifically, Figure 8A shows the raw StO 2 parameter map that was later preprocessed to enhance the contrast of blood vessels before using as an input to the B-COSFIRE algorithm for blood vessel segmentation. Figure 8B presents the binary blood vessel segmentation mask obtained by thresholding and postprocessing the B-COSFIRE response image and overlaid onto the StO 2 parameter map. A comparison between the applied binary mask in Figure 8B and RGB image in Figure 5 confirms that most blood vessels visible by the naked eye were accurately segmented, as well as some blood vessels hardly visible from the image. Blood oxygenation within the blood vessels was calculated by applying the binary mask and the skeleton mask to the StO 2 parameter image of subject 2. For 4T1 mammary carcinoma, the respective mean values of blood oxygenation in the tumor area and healthy tissue were 13.47% ± 18.51% and 25.52% ± 14.49%, as mentioned in Section 3.1. Notably, healthy tissue oxygenation is an average of the vessel oxygenation and the surrounding tissue oxygenation. With the blood segmentation masks available, we could separate the two areas. StO 2 values within the blood vessels and the healthy tissue excluding vessels were 37.49% ± 5.46% and 24.70% ± 14.55%, respectively. As expected, these results show that blood oxygenation within the vessels is much higher than in other areas of healthy tissue, and the blood vessel oxygenation variation is much lower.
Moreover, Figure 8C shows the blood vessel skeleton representing vascular architecture overlaid applied to the StO 2 map. From the skeleton map, the oxygenation profile within two annotated blood vessels was determined in the direction indicated by the red arrows and is presented in Figure 8D. Such an approach makes it possible to determine the changes in blood oxygenation along the vessels within the tumors or the surrounding tissue. Since the image pixel size is known, the length of the blood vessels could also be accurately calculated. For example, the F I G U R E 8 Blood vessel segmentation from the StO 2 map (subject 2) using the B-COSFIRE algorithm: A, StO 2 parameter map; B, postprocessed binary segmentation mask applied to StO 2 map; C, skeleton representing vascular architecture applied to StO 2 map and D, oxygenation profile along blood vessels 1 and 2 (annotated in B) length of the blood vessels 1 and 2 between the red markers was 0.44 and 0.65 cm, respectively. The oxygenation and blood volume in the blood vessels could also serve as tumor biomarkers.

| Tumor volume estimation
This section presents the tumor volume estimation based on three-dimensional optical profilometry integrated with the hyperspectral imaging system. First, tumors were segmented from the normalized hyperspectral images based on the SAM method Equation (2). The segmented tumor masks were then post-processed to fill the holes within the tumor and remove the outlying points to get the binary masks. Figure 9A shows the 525 nm band of the hyperspectral image of a TS/A murine mammary adenocarcinoma, whereas Figure 9B presents the binary segmentation mask of the tumor obtained with the SAM method. There is a good agreement between the visually estimated tumor location and its binary mask. Figure 9C shows the hyperspectral image overlaid onto the 3D point cloud extracted using the optical profilometry module, enabling 3D inspection of the imaged tumor.
The tumor volumes estimated from 3D OP data using the two approaches described in Section 2.6 are compared to the standard vernier caliper procedure in Table 1. All methods roughly agree; however, the volumes estimated from the profilometry data using the first approach (OP1 in Table 1) generally overestimate the manual measurements by 2.9% to 8.5%; on the other hand, tumor volumes estimated from the profilometry data using the second approach (OP2 in Table 1) are either smaller or larger by a minimum of 7.4% and a maximum of 26.3%, compared to the standard technique using a vernier caliper. The discrepancies between the approaches mainly result from the fact that the actual tumor shapes are not ellipsoids as assumed by the  standard approach. For example, a 4T1 mammary carcinoma presented in the RGB image ( Figure 5) has an irregular shape and is relatively shallow in the middle.

| DISCUSSION
This article aims to show that distinct characteristics of tumor tissue can be investigated and quantified using a multimodal optical imaging approach combining hyperspectral imaging and optical profilometry. The former has shown much promise in highlighting the physiological differences in tumors and healthy tissue, such as blood volume fraction and tissue oxygenation, and morphological differences reflected in the tissue scattering properties. Furthermore, blood vessels within the tumor and its environment can be segmented from the extracted blood parameter maps, allowing for monitoring angiogenesis in early-stage solid tumors and blood vessel growth over time. On the other hand, profilometry can be used to measure tumor size and thus monitor tumor growth by accurately estimating tumor volumes from three-dimensional point clouds. To the best of our knowledge, this is the first study to apply the combination of hyperspectral imaging and three-dimensional optical profilometry to examine murine tumor models in a preclinical pilot study. Figure 3 shows the normalized measured reflectance skin spectra of various tumor models and surrounding healthy tissues. The spectra allow direct discrimination between most tumor models (4T1, B16-F10, TS/A) and healthy tissue without needing to extract the tissue parameters from the hyperspectral images since the spectra change dramatically when the disease is present.
The hyperspectral images of murine tumor models were fitted using the IAD algorithm to extract tumor tissue's physiological and morphological parameters and the surrounding healthy tissue. Figure 4 shows the measured and fitted spectra of a 4T1 mammary colon adenocarcinoma (subject 2). Specifically, the spectra from within the less oxygenated, more oxygenated, and tumor tissue are presented. Here the characteristic shape of the spectra is altered in the tumor model, demonstrating that hyperspectral imaging is very sensitive to small changes in oxygenation, especially in the 550-600 nm spectral range. The fit results also show a good agreement between the measured and fitted spectra, suggesting that the IAD algorithm is accurate and robust in extracting the tissue parameters. However, our recent publication [57] showed that different sets of model parameters could yield precisely the same reflectance spectra, meaning there is already some intrinsic uncertainty in parameter estimation.
A comparison of the tissue parameters between all tumor models and the surrounding healthy tissues considered in this study is presented in box plots in Figure 7. It shows that physiological properties such as f m , f Hb , f HbO 2 , StO 2 and total hemoglobin content are generally higher in tumors than in healthy tissues. Meanwhile, the opposite applies to scattering properties such as a and b. Lower scattering parameters for tumors could be related to the structural changes within the diseased tissue, such as the local collagen fiber reorganization and alignment and cell and nuclei sizes [60]. However, the tumor oxygenation is expected to be lower than in the healthy tissue due to the increased demand for oxygen and disrupted supply, which opposes our findings [1][2][3]. A potential explanation is that the oxygen supply in small tumor models is sufficient to meet the growing tumor demand, as reported by Ziemer et al. [61]. On the other hand, large tumors are usually more hypoxic in the center and well-oxygenated in the tumor margins [62]. It is also possible that high oxygenation within the investigated tumors could be attributed to petechiae from skin irritation due to depilation. Generally, low oxygenation values are probably because the skin is typically among the least oxygenated tissues [63]. We have also found that total hemoglobin content is four times higher in tumors than in healthy tissue, meaning there is more blood within the diseased tissue due to an increased demand for oxygen and nutrients. Finally, we confirmed that the differences in tissue parameters between the tumors and healthy tissue are statistically significant (p < :05) for all parameters except a, as determined by the Mann-Whitney U-test.
We also investigated the oxygenation within the blood vessels segmented from the StO 2 map ( Figure 8A) using the B-COSFIRE algorithm. Figure 8B shows the binary segmentation mask overlaid onto the blood oxygenation map of a 4T1 tumor model, and Figure 8C shows the blood oxygenation of the vessel skeleton map. Using the B-COSFIRE, it is possible to isolate the oxygenation within the healthy tissue, excluding the blood vessels. We showed that blood oxygenation in the tumorfeeding blood vessels is approximately 50% higher than in the remaining healthy tissue and that the oxygenation variation is much lower in the vessels than in the tissue. Moreover, oxygenation profiles for two blood vessels were plotted in Figure 8D, allowing us to determine changes in oxygenation along the vessels. A combination of the presented information could be beneficial in monitoring the angiogenesis and blood vessel growth in tumors and the surrounding tissues as the disease gradually progresses and in monitoring the dynamic processes such as stroke. We also demonstrated that the vascular architecture and blood vessel length could be determined from the vessel skeletons, opening many possibilities to quantify vascular morphology, such as vessel skeleton density [64].
A standard approach to monitoring tumor progression in murine tumor models is the measurement of tumor volumes [46]. Usually, daily measurements of the three tumor axes are performed with a vernier caliper, and the tumor volume is calculated using Equation (7). However, approximating solid tumors with an ellipsoid yields errors in tumor volume estimation since the tumors are not ellipsoids. Our results show that the discrepancy between the manual approach and the estimated volumes from the profilometry data is within 26% ( Table 1) when tumor volume is estimated using Equation (8). However, when tumor volumes are estimated from profilometry data by measuring three perpendicular axes and using Equation (7), the differences with the manual measurements are less than 10% and as low as 2.9%. The latter approach is more precise because it uses the same method Equation (7) as the manual approach yet still overestimates the tumor volume. We believe this is because the vernier caliper slightly squeezes the tumor when performing manual measurements, resulting in lower tumor volumes.
On the other hand, the tumor volumes from the profilometry approach using Equation (8) are more different from manual measurements than the profilometry approach using Equation (7), presumably related to the actual tumor shapes differing from the ellipsoid shape. However, because of a remarkably high spatial resolution of the 3D OP module (0.1, 0.1 and 0.05 mm in the X, Y and Z directions), we firmly believe that tumor volumes estimated from profilometry data using the second approach Equation (8) are more accurate than two other measurements. We showed that with the help of the 3D OP, tumor volumes could be estimated more accurately than using a vernier caliper.
Furthermore, Figure 9C presents the spectral image overlaid onto the 3D point cloud. This visual representation demonstrates the power of combining the complementary information using different optical imaging modalities to improve the disease diagnosis by helping convey the relevant diagnostic information to the observer.
This pilot study demonstrated that a multimodal optical imaging system comprising hyperspectral imaging and three-dimensional optical profilometry could discriminate between tumor and healthy tissue based on quantitative information about tissue physiology, pathology and morphology. However, further studies of specific tumor models are needed with enough subjects (at least [10][11][12][13][14][15][16][17][18][19][20] to assess the tumor growth, angiogenesis and blood vessel growth characteristics with a high degree of statistical significance. A sufficient number of subjects for specific tumor models could also help identify that particular model's physiological and morphological characteristics and help discriminate different tumor types. These properties might include but are not limited to physiological parameters such as melanin volume fraction and tissue oxygenation. The former is usually high in melanoma tumors [65]; the latter is often low in tumors and even lower in highly aggressive phenotypes [3]. Different tumor models could also be discriminated based on morphological properties such as graininess of the surface or highly tortuous and porous blood vessel structure [3][4][5].
One limitation of this study is that, at the moment, only blood vessels in the tumor environment can be detected and segmented, especially in sizeable tumors. We believe this could be attributed to some parts of tumors being out of focus due to the shallow depth of field of the camera lens. Since the imaging plane was in the middle of the tumors in the Z direction, small blood vessels within the tumors could be blurred and thus not distinctive. Another possibility is that blood vessels in large tumors are located deep within the tissue and therefore not reached by the incident light in the 400-1000 nm spectral band.
Light penetration depth in the human skin in this spectral band is between 0.5 and 2.5 mm [66], but human skin is significantly thicker than murine skin. The human epidermis and dermis are generally 100-150 μm and 1-4 mm thick [6,67], respectively, compared to 10-50 μm and 0.1-0.7 mm in mice [67], meaning that light in the 400-1000 nm spectral band penetrates these two skin layers in mice and can even probe the subcutaneous adipose tissue layers. The absorption of the subcutaneous layer is defined by hemoglobin, lipids and water, but it was not considered in our two-layer model because absorption by water and lipids is especially relevant in the NIR spectral band for wavelengths higher than 1000 nm not used in this study. Moreover, Bjorgan and Randeberg [68] proved that by exploiting the scaleinvariance of the reflectance modeling, it is possible to consider just the epidermis and dermis without fully modeling complex deeper layers. Therefore, a two-layer murine skin model should sufficiently accurately model light propagation in biological tissue and estimate tissue parameters describing tissue physiology and morphology.
In our recently published research work [57], we showed that the accuracy and robustness of tissue parameter extraction by the IAD algorithm improved significantly as more model parameters were fixed. For example, the standard deviation of f Hb decreased eightfold as the number of free model parameters was reduced from 11 to 5. In the present study, all model parameters of the two-layer murine skin model were estimated.
Although we have no ground truth values for estimated tissue parameters, based on our previous observations for simulated reflectance skin spectra, the mean values of fitted parameters are within AE1 STD from the actual values. We also demonstrated that the correlation between the estimated parameters depends on the selection of fixed parameters [57]. In the present study, f m and d e exhibited the highest negative correlation (Pearson correlation coefficient was p 12 ¼ À0:7538), meaning a decrease in epidermal thickness yields a higher melanin volume fraction. It can be explained by the fact that melanocytes in the epidermis are the major producers of melanin, which is therefore the main absorber in the epidermis [6]. As the layer thickness decreases, the volume fraction of melanin must be increased to compensate for a change in light absorption. On the other hand, no correlation (p 12 ¼ À0:0708) was found for f m and f HbO 2 , meaning a change in one parameter does not affect the other at all.
A relevant limitation of the present study is low reflectance values for the B16-F10 melanoma tumor model ( Figure 3B) due to the high melanin volume fraction within the melanoma, f m ¼ 6:46% AE 2:68%, leading to the inaccurate extraction of the tissue parameters such as oxygenation. To increase the amount of signal, the current exposure time of 250 ms would have to be increased significantly. However, this could lead to overexposure of the surrounding tissue with much less absorption from melanin, so the exposure time should be carefully optimized to satisfy both cases.
An important issue is also an interplay between melanin and tumor necrosis. The melanin volume fraction is overestimated in necrotized tumors ( Figure 6A) because the spectra of melanin and necrotized tissue seem similar, and the necrosis optical data is not included in our current tumor model. Therefore, we plan to obtain histopathological samples of the necrotized tumor tissue and perform microscopic hyperspectral imaging to obtain the corresponding reflectance and transmittance spectra. This data will be used to determine the optical properties of the necrotic tissue. Then, the necrotized tissue could be incorporated into the tissue model to better discriminate between the actual necrosis and tissue with high melanin volume fraction.
Tissue properties such as f m , f Hb and f HbO 2 could be estimated from the hyperspectral images by calculating standard RGB channel intensities as in Figure 5 and determining their concentrations from these channels. Nishidate et al. [69] estimated tissue parameters based on the RGB images and found a good match between the proposed method and the actual results. For example, the estimated oxygenation values were within 10% of the actual, while the error was 3.46% for the melanin concentrations. Similarly, the group by Spigulis proved the concept of single snapshot multispectral imaging for skin chromophore mapping [70] and applied it to extract values f m , f Hb and f HbO 2 for skin malformations [71,72]. The main advantage of estimating tissue properties from RGB images is speed due to a time-consuming process of least-squares fitting of reflectance spectra to estimate properties directly from hyperspectral images. However, extraction based on the RGB channel intensities does not consider the whole spectrum of LED illumination but individual values, possibly missing some crucial information because various chromophores express characteristic features in different spectral bands. Also, this approach generally involves estimating the absorption properties of tissue but rarely the scattering properties, which could provide relevant information about the tissue scatterers, such as their size. However, in the future study, we will directly compare the tissue parameters extracted either from HSI or projected RGB images to assess the feasibility of using a more computationally demanding full HSI analysis.

| CONCLUSION
This preliminary study examined eight mice with four different subcutaneously grown tumor models employing combined hyperspectral imaging (HSI) and threedimensional (3D) optical profilometry (OP). We demonstrated that hyperspectral imaging could differentiate between tumor models and healthy tissue based on the shape of skin reflectance spectra and the physiological (melanin volume fraction, blood oxygenation) or morphological (scattering properties) tissue parameters extracted from the hyperspectral images. We also demonstrated the blood vessel segmentation using B-COSFIRE based on tissue oxygenation maps, which could be used to monitor oxygenation within the blood vessels. Moreover, vascular architecture could be determined using the skeletonization approach, and the oxygenation profiles along the blood vessels could be measured. Additionally, blood vessel length and volume could be easily determined from the binary segmentation masks or skeletons. Additionally, tumor volumes could be accurately determined from the 3D point clouds obtained with the OP module, eliminating the need for manual tumor measurements using a vernier caliper.
By exploiting the combined information about tissue parameters, blood vessels and tumor volumes, we provide a valuable tool to monitor tumor growth and disease progression in murine tumor models. Following the biological validation in murine tumor models, the presented imaging system could be translated into the clinical setting and applied to non-invasive and contactless tumor detection in human patients.

CONFLICT OF INTEREST
The authors declare no financial or commercial conflict of interest.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.