2.1.

Tissue Phantoms

The preparation of the SiliGlass phantoms was based on a report by Konugolu Venkata Seker et al.²⁰ The exact recipe for preparation, including rigorous sonification, degassing, and temperature-controlled curing of the phantom, is described in detail in a previous paper.²¹ The medium for the phantom was PlatSil SiliGlass (Polytek), and the dedicated black absorber was Polycraft Black Silicone Pigment (MB Fibreglass, United Kingdom). The scatterers were silica microspheres (No. 440345, Sigma-Aldrich) with an expected mean sphere diameter of 9 to $13\mu \mathrm{m}$, as specified by the manufacturer. We used two low-absorption phantoms (Aa and Da) and two high-absorption phantoms (Af and Df). Approximate values of reduced scattering and absorption coefficient values obtained from the original paper²⁰ are given in Table 1.

Table 1

Predicted approximate optical properties of the four manufactured SiliGlass tissue phantoms at 600 nm as given by the recipe.20

During the phantom manufacturing process, all the components were weighted, along with the residual material in the mixing containers and on mixing utensils, to obtain exact mass fractions of individual components. Here, we assume that the mass does not change during curing. In addition, the density of the spheres was measured by weighting and mixing with water of pure microsphere scatterer. The determined ratio of densities was ${\rho }_{\text{sphere}}=1.196(1\pm 0.025){\rho }_{\text{SiliGlass}}$. This ratio was used to convert between volume and mass fractions while accounting for the overall phantom density change due to the inclusion of scatterers. The dimensions of the produced phantoms before further processing conformed to the sizes of the urine cups (A5-50.20.18APL, Mikro + Polo, Slovenia) used for the production. The phantoms were thus cylindrical, with a diameter of 50 mm and heights of about 35 mm.

After manufacturing the phantoms, they were cut in half and wet-sanded for macroscopic imaging. This cut-and-sanding process prevented both imaging of edges of the phantom near the container, where the concentrations of constituent components could be off, and removed most of the surface texture artifacts due to cutting. The manual wet sanding was performed with a silicon carbide wet sanding paper (CP918A, VSM, Germany) with progressive grits P320, P500, and P1200. The sanding patterns were randomized to minimize sanding traces. After sanding, the phantoms were washed thoroughly and inspected for residue under a microscope. Half of the phantom was cut into arbitrary long square rods with a cross-section of $\sim 1\mathrm{cm}$ by 1 cm and stored for further microscopy imaging. The sanded phantoms were placed on the imaging table of the system with the sanded side facing the camera and imaged in the reflectance geometry.

For the microscopic imaging, thin slices of the phantoms were cut from the rods obtained as leftovers after phantom cutting for macroscopic imaging. The samples were slit using a fresh cardboard knife blade using a handheld microtome; during this process, the samples were supported in the microtome utilizing a piece of hard plastic to prevent bending in the direction of cutting and consequent crumbling. The slices were then placed on a clean microscope object glass and imaged in the transmittance geometry. The sample thickness was measured by focusing the microscope on the glass slide and setting the $z$-position reader to zero. The area that was to be imaged was focused, and the offset of the $z$-axis gave its precise thickness: $(770\pm 10)\mu \mathrm{m}$ for phantom Aa, $(580\pm 10)\mu \mathrm{m}$ for phantom Af, $(250\pm 50)\mu \mathrm{m}$ for phantom Da, and $(240\pm 60)\mu \mathrm{m}$ for phantom Df. Since the cut was imperfect, the difference in thickness between different sides of the sample was given as the uncertainty. The Da and Df samples were cut thinner than the Aa and Af samples, resulting in more significant thickness variations. The whole process of phantom preparation for imaging is schematically presented in Fig. 1.

Fig. 1

Workflow of phantom post-production. Phantoms were first cut along a line 1 cm from the bottom. The top part was used for imaging using a laboratory HSI system, while the bottom part was used for microscopy. The phantom for HSI was sanded with wet sanding paper on the cut side and imaged in reflectance mode. The phantoms for microscopy were first cut into square rods and then sliced into thin slices using a handheld microtome. Imaging was performed on a microscopy object glass in transmission geometry.

2.2.

Hyperspectral Imaging and Data Normalization

Hyperspectral imaging (HSI) was performed on two size scales to examine the effects of the phantom microstructure on its optical properties, described in detail below.

Average reflectance spectra of tissue phantoms were obtained using our custom-developed macroscopic laboratory HSI system.²⁴ In brief, the system is based around an imaging spectrograph (ImSpector V10e, Specim, Finland) and a custom broadband LED light source enabling operation in the wavelength range between 400 and 1000 nm with a spectral resolution of 2.9 nm. Multiple different objectives can be used with the system, which results in the spatial resolution of the system at 500 nm of 0.3 and 0.1 mm for 17 and 50 mm lenses, respectively. The system employs a push-broom methodology, whereas the reflectance signal is acquired along a line scanned perpendicularly over the sample. The measured intensity is divided by the reflectance recorded on a white reference standard (WS1, Labsphere Inc.) with a 99.9% reflectance, obtaining normalized reflectance $R_{x,y}(\lambda )$ in the form of

The phantoms were also examined on a micro-scale using thin phantom slices where only a handful of scattering events occurred. To this end, our custom hyperspectral microscope was used.²⁵ The microscope operates in the filtered regime, acquiring a sample transmittance at one wavelength at a time. Scanning the whole wavelength range between 450 and 750 nm enables the acquisition of an entire hyperspectral cube with a spectral resolution of 2.5 nm and a maximal attainable spatial resolution of $1.3\mu \mathrm{m}$ with a 50× magnification objective. The system uses a custom-developed monochromatic light source based on a Czerny–Turner configuration monochromator with built-in on-line spectroscopic monitoring. All images were recorded with a spectral step of 1 nm and a camera integration time of 1 s; the images were binned in each spatial direction by a factor of 2. After the binning, the detector’s dark current was subtracted from the sample and reference measurements. The sample measurements were divided by the unobstructed beam reference image acquired using the same system settings, resulting in the normalized transmittance hyperspectral microscope images:

2.3.

Microsphere Size Distribution Determination

Microscopic properties of the scattering elements in the phantoms govern the overall scattering properties. The main ingredients of the phantom that contribute to the scattering are the glass microspheres. While the pigment particles also scatter light, their concentration is typically much smaller than that of the glass spheres. In addition, air bubbles that remained trapped in the polymer could potentially contribute to scattering, but their presence was minimized as best as possible through rigorous vacuuming of the samples.

Microscopic properties of the spheres were assessed by a scanning electron microscope (SEM)²⁶ Helios NanoLab 650 (FEI Company). A small amount (tip of a spatula) of glass sphere powder was placed on an adhesive imaging substrate and imaged at different magnifications; raw SEM images were preprocessed by locally enhancing the contrast using a MATLAB function adapthisteq. Afterward, circular objects within the image were detected using a MATLAB function imfindcircles with phase detection method, and their diameters were extracted. The spheres’ diameter was converted from pixels to $\mu \mathrm{m}$ using the scale bar included in the images.

2.4.

Absorption Properties Measurement and Determination

Absorption properties of both the polymer material and absorption component were measured in the transmittance mode by placing the material in a standard transparent polystyrene cuvette (Makro PS cuvette, 2711110, Ratiolab, GmbH) on a laboratory spectrometer (Lambda 960, Perkin-Elmer). The absorption coefficient was calculated from the transmittance values following the procedure described by Li et al.²⁷ The approach accounts for reflections from the cuvette walls by accounting for the medium and cuvette refractive index and absorption coefficient of the cuvette material. The refractive indices for polystyrene and SiliGlass were obtained from the literature.^21,28,29 The absorption coefficient of the polystyrene was measured by using the transmittance model on an empty cuvette with the expected content transmittance of 1.

2.5.

Calculation of Scattering Coefficient

The scattering coefficient and anisotropy were calculated for full and hollow quartz glass spheres using Mie theory.³⁰ A C programming language routine by Prahl³¹ was used for full spheres, and a routine from Bohren and Huffman in Fortran programming language BHcoat³⁰ was used for hollow glass spheres. Hollow glass spheres were simulated as a quartz spherical shell surrounding an air bubble. For all the simulations, the refractive index of the surrounding material was set to that of SiliGlass, as obtained from the literature.²¹ The refractive index of air was taken to be $n=1$, and the refractive index of ${\mathrm{SiO}}_2$ was obtained from the literature.^32,29 The result of the simulation is the scattering efficiency $Q_s$ of a sphere, which is related to the scattering cross-section of a particle $C_s$ through its geometrical cross-section $S$ as

In addition to the scattering coefficient ${\mu }_s$, the scattering anisotropy $g$ was also calculated. The BHcoat routine used for the hollow spheres does not support the calculation of the scattering anisotropy, so it was approximated by simulating just the air bubble (the hollow center of the sphere) surrounded by the SiliGlass, assuming similarity of refractive indices between SiliGlass and silica.

2.6.

Inverse Problem Solver

With known optical properties of pure phantom constituents, extracting volume fractions of individual components from the transmittance or reflectance values is possible by solving an inverse problem. To this end, a modified version of the Monte Carlo multi-layer (MCML) on CUDA^33,34 was used. A set of initial values was assumed, and a nonlinear fitting routine lsqnonlin in MATLAB (Mathworks) was used to compute the volume fractions of individual components iteratively. The reflectance values were fitted at 166 wavelength points between 420 and 980 nm. To evaluate the goodness-of-fit of the final results, three metrics were calculated based on the measured reflectance spectrum $R^{\text{measured}}({\lambda }_i)$ and Monte Carlo reflectance spectrum $R^{\mathrm{MCML}}({\lambda }_i)$. First, root-mean-square-error

2.7.

Software Environment

The simulations, including routines for calculating the scattering coefficient and inverse Monte Carlo, were performed on a desktop computer with an Intel i5 CPU, 16 GB of RAM, and Nvidia GTX1050-Ti graphics cards under the Arch Linux operating system. The C and Fortran code was compiled using the GCC 10.2.0 compiler collection. CUDA implementation of Monte Carlo was compiled for CUDA 11.1, running under the driver version 455.45.01. The preparation of the simulation input files and parsing of the simulation outputs was conducted using Python 3.9.1, which was running a web server implemented in Flask that enabled the invocation of individual computation routines over HTTP requests. The data processing was performed on a laptop computer (Intel i5, 16 GB of RAM) using MATLAB R2020a (Mathworks) with Optimization Toolbox and Image Processing Toolbox.

3.1.

Absorption Properties

The absorption properties were measured for the 1:2272 diluted pigment used to prepare the phantoms (termed the absorbing component or pigment in the continuation) and polymerized SiliGlass phantom medium. The absorption coefficients for both the pure polymerized polymer and pigment calculated from transmittances measured on the laboratory spectrometer are shown in Fig. 2. The values for absorption were compared to values given by Konugolu Venkata Seker et al.²⁰ Series $a$ phantoms have a reported ${\mu }_a(600\mathrm{nm})\approx 0.1{\mathrm{cm}}^{-1}$. For our measurements, mass fractions of diluted pigment were about $(3.7\pm 0.2)\%$ resulting in ${\mu }_a=(0.11\pm 0.01){\mathrm{cm}}^{-1}$, assuming the densities of the whole phantom are close to the absorber component. For phantoms $f$, the literature reports ${\mu }_a(600\mathrm{nm})\approx 0.9{\mathrm{cm}}^{-1}$. In our measurement, mass fractions of the diluted pigment were $(33\pm 1)\%$ resulting in ${\mu }_a=(0.99\pm 0.03){\mathrm{cm}}^{-1}$. From this comparison, the measured absorption coefficient agrees well with the previously reported approximate numbers within the uncertainties.

Fig. 2

Measured absorption coefficients of (a) SiliGlass polymer and (b) diluted pigment.

3.2.

Microsphere Size Distribution

The glass microspheres were imaged under SEM at four different magnifications, 80× to 3000×; in total, 12 images of different regions were acquired for microsphere size analysis. As expected from the manufacturers’ specifications, the images reveal a distribution of spheres of various sizes with some larger aggregates. Diameters of 8156 microspheres were determined from SEM images; an example of the process of diameter determination is shown in Fig. 3(a). To account for a large diameter variability among microspheres, images at different magnifications were used. Some spheres remained undetected, but visual inspection revealed that the detection probability was not related to the sphere’s size but to the sphere’s apparent brightness in the image. Based on this observation, undetected spheres should not introduce bias into the size distributions.

Fig. 3

(a) An example of microsphere diameter detection from SEM images (400× magnification); (b) microsphere diameter distribution normalized both in terms of number and volume fraction.

The microsphere size histograms were normalized to unity to obtain number fraction distribution. The volume fraction distribution was calculated by multiplying the number fractions by appropriate sphere sizes. The distribution was then normalized to unity and smoothed using a moving average filter in MATLAB, removing statistical noise due to the relatively small sample size. The resulting number and volume distributions are shown in Fig. 3(b), exhibiting two distinct peaks. Most spheres were relatively small, while a smaller number of larger spheres contributed prominently to the volume fraction. The average diameter of microspheres was $(13.47\pm 5.98)\mu \mathrm{m}$ when considering volume fraction distribution, which was close to the range of 9 to $13\mu \mathrm{m}$ given by the manufacturer.

The spheres could be hollow, so the wall thickness was introduced as another vital parameter besides the diameter. The wall thickness was measured by cutting four spheres of different diameters (3.4, 8.4, 18, and $24.6\mu \mathrm{m}$) in half using a focused ion beam (FIB)³⁵ integrated into the SEM device. After cutting the spheres, the resulting cross-section was imaged using SEM, with an example shown in Fig. 4. The smaller sphere ($3.4\mu \mathrm{m}$) was observed to be solid, while the largest sphere ($24.6\mu \mathrm{m}$) had a small air bubble that was not centered. For the medium-sized spheres (diameters 8.4 and $18\mu \mathrm{m}$), an average wall thickness $(0.9\pm 0.2)\mu \mathrm{m}$ was observed. The number of spheres cut in half was limited due to the long time needed to cut the sphere in half and the limited time on the device.

Fig. 4

SEM image of a sphere cut in half by an FIB. A microsphere with a diameter of $18\mu \mathrm{m}$, cut in half, is shown in the top and side view; magnification is 4000×.

3.3.

Predicted Scattering Coefficient and Anisotropy

Based on the determined scatterer microstructure, Mie theory numerical calculations were performed. First, some representative examples from the distribution were selected and are shown in Fig. 5. In general, smaller spheres exhibited more pronounced changes in scattering coefficient across the observed spectral range, as expected from Mie theory. The difference between the hollow and full spheres was particularly interesting for larger spheres. Both full and hollow spheres exhibited similar scattering coefficient values with more significant deviations from the average value for hollow spheres.

Fig. 5

Scattering coefficient and anisotropy for some example spheres calculated using Mie theory: (a) scattering anisotropy of small solid glass spheres, hollow spheres consisting of an air bubble with a quartz glass shell, and air bubbles (e.g., hollow spheres without the glass shell) submersed in the SiliGlass polymer; (b) scattering anisotropy calculated for full glass spheres and air bubbles in a SiliGlass polymer. The size value in the legend is the radius of the simulated sphere. In the cases of hollow spheres, the simulated shell thickness is $1\mu \mathrm{m}$ so that the air volume is the same as in the case of the air bubble without the shell.

A model that approximates the realistic properties of the scattering component was developed based on the measurements of sphere wall thickness. First, the wall thickness distribution was approximated to be box-shaped between 0.7 and $1.1\mu \mathrm{m}$ with a step size of $0.1\mu \mathrm{m}$. For each wall thickness, scattering coefficients and anisotropies were calculated following the size distribution of spheres. For the calculations, the distribution was interpolated to $0.1\mu \mathrm{m}$ step size and renormalized to unity. As a verification, a denser distribution was also examined but showed no significant changes in the final scattering coefficient. The spheres with a radius smaller than the wall thickness were simulated as full glass spheres (denoted in the equations by omitting the wall thickness parameter $d_{\text{wall}}$). Due to the limitations of the numerical codes, scattering anisotropy was calculated for full glass spheres in the case of small spheres and as air bubbles with radius decreased by the wall thickness in the case of larger spheres. The whole process is summed up by the model equations:

An approximation of both scattering coefficient and scattering anisotropy for the size-distributed silica microsphere scattering component was calculated using the model equations and is shown in Fig. 6 for 10% volume fraction of the spheres.

Fig. 6

(a) Scattering coefficient and (b) scattering anisotropy calculated for the determined size distribution of spheres in a non-interfering regime, accounting for the distributed wall thickness. The values are calculated for a 10% volume fraction of the spheres.

The scattering coefficient values can be again related to the values reported in the literature.²⁰ Given the 10% volume fraction of the scattering component in the simulations, the scattering coefficient for 1% sphere volume fraction would be ${\mu }_s(600\mathrm{nm},1\%)=31{\mathrm{cm}}^{-1}$ with the scattering anisotropy of $g=0.824$. The mass fractions of the scattering component determined during the preparation were $(1.7\pm 0.1)\%$ and $(7.5\pm 0.1)\%$ for phantoms in series $Ax$ and $Dx$, respectively. Correcting for the densities of the spheres, we arrive at volume fractions of $(1.5\pm 0.1)\%$ and $(7.3\pm 0.1)\%$ for both phantoms. Estimating the reduced scattering coefficient as ${\mu }_s^{\prime }=(1-g){\mu }_s$ gives values ${\mu }_s^{\prime }(600\mathrm{nm})=8.2(1\pm 0.1){\mathrm{cm}}^{-1}$ and ${\mu }_s^{\prime }(600\mathrm{nm})=39.8(1\pm 0.1){\mathrm{cm}}^{-1}$ for phantoms in $A$ and $D$ series, respectively. These results agree with the ones reported in the literature to the order of magnitude but exhibit deviations that could be due to the specifics of the hollow particle scattering phase function. For example, changing the scattering anisotropy to $g=0.9$ would cause the reduced scattering values to agree with those from the literature given in Table 1 perfectly.

Furthermore, seemingly exaggerated spectral features were observable between 400 and 500 nm. These could be attributed to a rather similar glass shell thickness in the case of hollow spheres. In general, if the shell thickness variations are larger than accounted for in this study, this would represent about 7% relative uncertainty in the presented data.

3.4.

Determination of Component Volume Fractions

The knowledge of sphere microstructure enabled numerical calculation of the predicted scattering coefficient. The measured absorption properties of the phantom constituents made it possible to determine the component volume fractions based on the HSI-measured spectra. The four phantoms were measured using the HSI system. Spectra were averaged from 20 points in the center of the phantom, as shown in Fig. 7(a). Average reflectance spectra, normalized to a white reference, are shown in Fig. 7(b) as colored lines with noticeable oscillations around 500 nm and between 800 and 900 nm. A comparison of reflectance spectra with the numerically calculated scattering coefficient in Fig. 6 confirms that similar shapes are observed.

Fig. 7

Measured and fitted spectra of tissue phantoms using a laboratory HSI system. (a) RGB projections of phantoms with sampling points used to calculate average reflectance spectra are denoted. Please note that the apparent inhomogeneity in the images results from RGB projection and aggressive histogram equalization in the visualization software. (b) The calculated average reflectance spectra for all four phantoms are displayed (colored lines) along with the results of inverse Monte Carlo fitting (black lines). Standard deviations of spectra for the sampled points are smaller than the line thickness.

Beyond qualitative comparison, the known sample microstructure made it possible to determine absolute volume fractions of phantom components. This was achieved by fitting the measured average reflectance spectra using the MCML on CUDA algorithm. The measured absorption and calculated scattering spectra were used as inputs, while the volume fractions of constituent components were left as free parameters. Graphically, the results of this inverse problem are shown in Fig. 7(b) with full black lines, achieving a good agreement between the measured and fitted spectra. Note that the oscillations around 500 nm are slightly exaggerated in the fit results, which is attributed to the limited knowledge of sphere wall thickness. Numerically, the results also exhibit good agreement between the experimental and theoretical values, as indicated by the deviations of the fit from the measured values and low RMSE numbers in Table 2. Numerical results for volume fractions of scattering and absorption components obtained by this process are presented in Table 3 and again show good agreement within the uncertainties of the method.

Table 2

Evaluation of goodness-of-fit for Monte Carlo fitting. For each phantom, the absolute mean deviation of the reflectance, as well as the maximal absolute value deviation, is given for two wavelength windows. In the last column, the root-mean-square-error of the fit for the whole range between 420 and 980 nm is reported.

Table 3

Values of mass fractions obtained by weighting during the preparation for absorber wabs and scatterer wscat in comparison to the mass fractions of the absorber w˜abs and scatterer w˜scat, as obtained by solving the inverse problem using computed scattering and measured absorption coefficients with the Monte Carlo method.

3.5.

HMI Evaluation of Tissue Phantoms

First, RGB projections of the normalized images were calculated by projecting the transmittance spectra to CIE XYZ color space and applying a standard D65 illuminant; the resulting images are shown in Fig. 8(a). The larger spheres visible in the images were dispersed equally over the sample, but aggregates of spheres, which were already observed in SEM images, are also visible using the HMI as irregularly shaped small inclusions. It appears that despite the rigorous sonification and mixing during the phantom preparation, not all the clusters of microspheres were successfully broken.

Fig. 8

Microscopic HSI of phantoms. In column (aX), RGB projections of the hyperspectral image (10x magnification) are shown. In column (bX), the spectra sampled at the sampled points denoted in column (aX) are shown. The first row (a1), (b1) shows phantom Aa, the second row (a2), (b2) phantom Af, the third row (a3), (b3) phantom Da, and the fourth row (a4), (b4) phantom Df. Locations 1 to 3 on each image correspond to the dark areas with absorber agglomerates, while areas 4 to 7 are sampled on the brighter spots near the scattering spheres. Spectral features due to the Mie scattering are observable in feature-rich spectra in those points. All spectral plots share the same wavelength axis at the graphs’ bottom. All spectra are binned over 16 pixels of the original image centered around the point marked in (aX).

The transmittance spectra were sampled from the hyperspectral images at multiple sites to analyze different phantom parts [Fig. 8(b)]. The RGB images and corresponding spectra revealed two distinct area types of the phantom. The locations of darker and brighter areas are denoted in Fig. 8 with 1-3 and 4-7, respectively.

The darker areas, rich with the absorbing material, exhibited low transmittance, and it is evident that the ink particles were not dispersed equally over the whole phantom but were agglomerated. Only minor differences in transmittance values between high-absorption and low-absorption phantoms were observed in these darker areas; furthermore, transmittance spectra were practically constant in these areas. Both observations indicate that most of the light in these areas was absorbed over the whole wavelength range. It is important to note that under close examination, the darker areas corresponded to the rougher parts of the sample in some cases. This raises two important conclusions: (1) the presence of the pigment introduced a defect in the silicon matrix that causes the phantoms to break more frequently along areas rich in the absorber, and (2) the results presented for the darker areas could be, to varying degrees, also impacted by the uneven surface.

In brighter areas, the transmission spectra for each phantom exhibited a strong dependence on the wavelength, indicating the contribution of single scatterers to the spectrum. Since the thickness of the sample was small, the number of scattering events was likely low. Therefore, the characteristic Mie shape of scattering from a single particle instead of the average scattering was to be expected, which was confirmed by the oscillatory spectra.