1. Introduction

Brain tumors are divided into benign and malignant tumor entities [1]. Regardless of the classification, surgical removal is widely considered the preferred treatment approach. Complete resection leads to prolonged patient survival, especially in the tumor group of malignant tumors such as glioblastoma and metastases [2–5]. A significant challenge in surgery is to differentiate intraoperative between healthy and diseased tissue in the central nervous system and to perform tumor resection without damaging surrounding essential brain areas. Numerous pre- and intraoperative neurosurgical imaging techniques are commercially available to assist the surgeon during surgery. Some of the most used intraoperative modalities are neuronavigation, fluorescence-guided tumor resection, intraoperative magnetic resonance imaging, and intraoperative ultrasound [6–8]. All the systems listed are error-prone and disrupt the intraoperative workflow, mainly when using magnetic resonance imaging and intraoperative ultrasound. Neuronavigation depends on preoperative images, which may turn inaccurate if anatomical changes occur during surgical tumor resection (brain shift). Another option for experienced surgeons to distinguish between healthy and diseased tissue is to use intraoperative haptic information transmitted by microsurgical instruments. This haptic information reproduces the elasticity behavior of the tissue. Elasticity is of great clinical interest since tumor tissue often behaves mechanically differently than healthy tissue [9]. For example, during breast self-examination, the sensitivity of the hands, which is approximately 40 µm at the fingertips [10], is used to detect breast abnormalities. However, it is a subjective technique that depends highly on the person performing the procedure. For this purpose, various automatic elastography techniques have been used and tested over the years in neurosurgical settings to map the elastic behavior of the tissue, most notably intraoperative ultrasound elastography and magnetic resonance elastography in a preoperative setting [11–13].

An alternative imaging technique is optical coherence elastography (OCE). It is based on the principle of optical coherence tomography (OCT), an interferometric measurement technique for non-invasive acquisition of volumes and cross-sectional images [14–16]. Swept-source OCT (SS-OCT) systems are based on a narrow-band, tunable light source whose high sweep rate and wide spectral bandwidth enable remarkably high imaging speed and axial resolution [17]. Following the introduction of the first Fourier Domain Mode Locked (FDML) laser, continuous improvements have enabled SS-OCT systems to operate at A-scan rates of up to 5.2 MHz with single-spot and 20 MHz with parallel four-spot imaging [18,19]. For most OCE techniques, resolution ranges from 15 to 100 µm with imaging depths between 0.5 and 3 mm [14]. However, there are even higher-resolution systems. Quince et al. [20] showed that even a resolution of approximately 3 µm could be achieved with OCE. The basic concept of OCE is to evaluate the response of a sample to a mechanical load [14,16,21–23]. Standard methods include compression, acoustic radiation, and air or laser pulses [14–16]. For instance, “optical palpation” methods are being researched as a possible alternative to subjective manual palpation in breast cancer diagnostics [24,25] or for measuring soft contact lenses [26]. Subsequently, the corresponding sample response is usually analyzed by speckle-tracking techniques [16,21,23] or phase-sensitive approaches [16,22,27].

Speckles, or random sub-pixel intensity fluctuations in the OCT image, occur when coherent light is back-reflected from rough surfaces. It is traditionally used in OCE to analyze the mechanical properties of the tissue by tracking the speckle pattern of a predefined multi-pixel window using image correlation. However, to allow such correlation, the speckle field must remain comparable before and during displacement. Therefore, the maximal possible sample movement is limited between two frames and must stay within a certain threshold [15,28]. At the same time, a sufficient mechanical load must be applied to the sample to allow speckle motion analysis. But because of its dependence on the speckle size, the resolution is lower than with phase-sensitive OCT (PhS-OCT). While speckle-tracking determines intensity changes in the OCT image, PhS-OCT analyzes the complex phase signal obtained by Fourier transforming the detected interference fringes. Since the resolution depends on the applied wavelength and OCT relies on light sources with wavelengths of approximately one micrometer, stiffness variations in the sub-micrometer range can be determined by measuring the phase shift of the backscattered light. This results in a high displacement sensitivity, one of the main advantages.

A limiting factor arises because the phase is restricted to $2\pi $. Therefore, the maximum measurable displacement between two A-scans can only range between $- \pi $ and $\pi $. This leads to an ambiguous phase information, also called wrapped phase. Those $2\pi $ ambiguities must be identified and corrected during a phase unwrapping process. Due to the delicate nature of the brain tissue and the cerebrospinal fluid surrounding it, the maximum load capacity is inherently limited. In addition, brain tissue's nonlinear and viscoelastic properties restrict the applicability of conventional mechanical calculations to the linear-elastic part of the stress-strain curve, necessitating modest forces. Visualization 1 illustrates the intricate behavior that occurs when excessive force is applied. The force is strong enough to cause a vigorous shaking motion, challenging elasticity determination. Therefore, there is a great need to measure smaller displacements, reinforcing the need for PhS-OCT. However, a sufficiently high load must be applied to the sample to obtain an accurate displacement curve and allow further mechanical analysis. Moreover, for in vivo brain measurements, a significant load is necessary to facilitate subsequent corrections for patient movements related to heartbeat or respiration. Since we are working with a MHz OCT system, analyzing the phase between successive A-scans is not advisable. An imaging speed reduction would be necessary as the phase shift is insignificant. However, the issue of phase wraps becomes particularly crucial when considering inter-frame phase analysis due to the extended time intervals between two referenced voxels. In addition, brain tissue is one of the most intricate and elastic tissues in the human body. They become inhomogeneous during tumor development and exhibit considerable differences in composition, structure, and mechanical properties. These differences pose a challenge for the unwrapping process. Despite the introduction of numerous unwrapping algorithms, these often have problems correcting rapid, large inter-frame variations.

To overcome these challenges, this study presents and evaluates a one-dimensional unwrapping approach. Measurements were done using a MHz OCT imaging system and a self-built non-contact air-jet device [29]. Building on a previously introduced method [30,31], the process was further developed and enhanced to improve its robustness and time efficiency. The detection and correction of multiple 2π wraps between successive B-scans occur within an approximate processing duration of 0.11 s per frame on the GPU. This method enables precise and reproducible sample displacement information for each pixel in depth and over the entire acquisition time. The algorithm analyzes the signal derivatives and initiates further correction and evaluation steps based on the outcome. These steps include applying an environment-aware predictive approach to examine the impact of the estimated phase value on the phase progression. Final local outliers are detected and rectified during post-processing using a sliding 5 × 5 median window. Within our datasets, the algorithm demonstrates the capacity to assess inter-frame phase shifts of up to 4π and an overall sample displacement of approximately 18 µm. This wide range of detectable shifts allows for comprehensive visualization. The results are then displayed through two-dimensional maps that visualize the tissue displacement for each pixel at a specific point in time. We extensively verified and validated the algorithm's reliability and performance using a series of phantoms. In addition, the results of ex vivo studies on human brain tumors were analyzed and supported by histological analyses.

2. OCE imaging system

A commercially available SS-OCT imaging system was used to perform OCE measurements (OMES, Optores GmbH, Germany). The system was equipped with a 3.2 MHz FDML laser operating at a central wavelength of 1310 nm, a bandwidth of 100 nm, and a measured laser power of 17 mW on the sample. The scan head was equipped with a scanning objective lens with an effective focal length of 54 mm, a numerical aperture of 0.037, and a working distance of 42.3 mm (LSM04, Thorlabs Inc., USA). The system has an axial resolution of 15 µm and a lateral resolution of 21 µm. An analog-to-digital converter with 4 GSample/s and a 12-bit sampling resolution was used for data acquisition (ATS9373, AlazarTech, Canada). All data processing was performed on a GPU-based High-Performance-Computing (HPC) server system equipped with two AMD CPUs (AMD EPYC, Advanced Micro Devices Inc., USA) and an Nvidia graphic card (Nvidia Quadro RTX A5000, NVIDIA Corporation, USA). The measurements were performed on a vibration-isolated optical table. To investigate the temporal phase response, 4 mm line scans were conducted with unidirectional scanning using the fast X-galvanometer scanner. In total, 792 frames with 512 A-scans were recorded. This corresponds to a B-scan rate of 2.45 kHz, equivalent to a measurement duration of 0.32 s.

As shown in Fig. 1, a self-developed air-jet system (Medizinisches Laserzentrum Lübeck GmbH, Germany) was used for initial experiments, which can be operated at adjustable flow rates and emits forces between 1 µN and 40 mN [29]. Two 1.4 mm nozzles were attached to the scan head at an angle of 45° to the surface and a sample distance of approximately 10 mm. During the assembly process, particular attention was taken to ensure that the peak of the Gaussian profile was in the center of the B-scans. The Gaussian force profile in air has a full width at half maximum (FWHM) of 2 mm. A complete measurement examines the effect of a single 200 ms air pulse. Since we use such short pulses, the air-jet does not cause any dehydration of the sample. For PVA phantoms, an air-jet force of around 30 µN was applied, while for the tumor measurements, the force varies between approximately 30 and 85 µN, depending on the properties of the tissue. To identify suitable force values for different samples, it is important to progressively increase the force until the entire displacement curve, containing all features, can be unwrapped. Once the optimum force value has been determined, it usually applies reliably to all samples of the same type.

 

Fig. 1. Experimental setup for tissue indentation and elastography measurement. Left: Scan head with two 1.4 mm nozzles attached to the objective from which the air pulses are ejected. Right: Air-jet system with exemplary force curve.

Download Full Size | PDF

2.1 Samples

Initial tests and validations of the unwrapping algorithm were performed on two different sample types, including tissue-mimicking phantoms made of polyvinyl alcohol (PVA) (Söring GmbH, Germany). Such PVA phantoms were subjected to repeated freezing and thawing cycles. Each cycle increased the branching of the polymer chains, which in turn increased the mechanical stiffness of the phantom. An experienced neurosurgeon evaluated a series of homogeneous phantoms to assess their haptic properties and elasticity by manual palpation. This process led to the development of a series of phantoms encompassing healthy and tumorous brain tissue properties, ranging from 875 Pa for healthy tissue to 16.38 kPa for tumorous tissue. Our project partner measured the Young’s modulus of the phantoms within the research project using unconfined compression. Also, phantoms with stiffness inhomogeneities were produced to simulate the transition between healthy and tumor tissue. Here, more rigid PVA inclusions were embedded in a softer PVA surrounding. Their resemblance to healthy and tumorous tissue was verified through palpation by an experienced neurosurgeon.

In addition, studies on ex vivo meningioma, glioblastoma, and metastatic human brain tumor samples were performed (Study No.: AZ 19-319, approved by the Ethics Committee, University Hospital Schleswig Holstein, Campus Lübeck). After surgical removal, the tumor was placed directly into a sample container filled with Ringer’s solution consisting of 8.6 g NaCl, 0.3 g KCl, and 0.33 g CaCl₂ per 1000 ml distilled water. It was then placed in an ice-filled cooler box. 3D-printed molds containing two cylindrical elevations with a diameter of 12 mm were used to ensure a uniform shape of the tumor samples and to allow accurate correlation with histologic evaluations. These molds were placed in biopsy embedding cassettes (Sanowa Laborprodukte GmbH, Germany) and filled with liquid agar. After the agar had cooled, the molds were removed. The tumor samples were cut into 12 mm segments with a punch and placed in these wells for future correlation of the OCE data and the histological sections. After imaging, the left half of the acquired line scan was labeled with green histological dye (WAK HM G 1/60, WAK Chemie Medical GmbH, Germany). Subsequently, the samples were embedded in paraffin and stained with hematoxylin & eosin (H&E). Ten sections were taken along the marked line at 200 µm distances to ensure a histological image of the acquired line scan and to enable histopathological evaluation. The slides were then digitized (VENTANA iScan HT, Roche Holding AG, Switzerland), registered, and correlated with the acquired OCT images. A detailed description of the registration process can be found in the study by Strenge et al. [32]. Since the green histological dye is only visible on the left side of the sample surface, the orientation of the histological image could be determined and retrospectively correlated with the sample position during acquisition.

3. Investigation of the phase signal

As any sample motion impacts the phase signal in the OCT interference pattern, it can provide high-resolution information about the mechanical response of the sample. By performing a Fourier transform, a complex signal can be generated, resulting in a depth-resolved reflection profile (magnitude of the complex signal) and the phase angle ${\phi _{x,y}}(t )\in {\mathrm{\mathbb{R}}^{M \times N \times L}}$ of the interference signal. The obtained unwrapped phase ${\varphi _{x,y}}(t )\in {\mathrm{\mathbb{R}}^{M \times N \times L}}$ for each frame t and each pixel $({x,y} ) $ can then be expressed as:

However, in reality, many factors influence negatively the real achievable value. Therefore, to investigate the phase noise ${\phi _{noise}}$ of our experimental 3.2 MHz OCT setup, a measurement was performed on a stationary PVA phantom without applying any mechanical load. Subsequently, the phase data was unwrapped using the standard MATLAB unwrap function (The MathWorks, Natick, Massachusetts, USA). The results are presented in Fig. 2. The left part of Fig. 2 shows the phase progression of different pixel positions within the central A-scan, spaced 40 pixels apart. The blue curve corresponds to the lowest depth, while the green curve represents the highest depth. The right part of Fig. 2 illustrates the phase difference, denoted as $\Delta {\phi _{noise}}$, between the first and all subsequent depths. It simulates the situation where a reference reflection is used to unwrap the phase information. The known phase change caused by the reflection can be used as a reference or calibration value to unwrap the signal, as already demonstrated by other groups. The phase noise remains similar in the first three depths but becomes higher with increasing depth of the examined phase signal. This trend can also be observed for $\Delta {\phi _{noise}}$. This is because the signal intensity decreases with increasing sample depth. Nevertheless, the standard deviation remains at a comparatively low value of 10 and 30 nm. It should be noted that the phase noise of the entire system is considered. Thus, all influencing factors like environmental movements, trigger jitter, electronic noise, and clocking instability are included. Many additional factors influence the real phase instability, making the phase noise value higher than the mathematically calculated ideal value. However, a high phase stability of the system is an essential requirement for obtaining displacement maps, regardless of the correction algorithm applied. The results indicate that the MHz OCT setup has a remarkably low $\Delta {\phi _{noise}}$ value, indicating no global timing jitter, which is usually a problem with SS-OCT systems.

 

Fig. 2. Analysis of the phase instability of the MHz OCT setup. One measurement was taken on a PVA phantom without applying any mechanical load. The left part shows the phase progression of different pixel positions within the central A-scan, spaced 40 pixels apart. The blue curve corresponds to the lowest depth, while the green curve represents the highest depth. The phase difference between the first depth and all subsequent ones is shown on the right side. For better visualization, artificial offsets were set. STD: Standard deviation.

Download Full Size | PDF

3.1 Phase unwrapping algorithm

This chapter introduces a one-dimensional unwrapping algorithm that works analogously to the interpretation of a displacement-time graph. It recognizes abrupt changes in velocity and acceleration over time, allowing effective phase signal processing. A technique is then applied that corrects the phase information using an environment-aware predictive approach to assess the impact of the corrective actions on the future phase values. The phase between two referenced A-scans from subsequent B-scans is evaluated throughout these processes. Figure 3 gives an overview of the phase signal processing steps. Code 1 [37] contains the MATLAB R2023a implementation of the algorithm with a processing time of approximately 0.11 s per frame on the GPU.

 

Fig. 3. Flowchart of the phase unwrapping algorithm. To enhance comprehension, the plots on the right depict various exemplary datasets, corresponding to the steps highlighted in green in the flowchart.

Download Full Size | PDF

After the phase data is extracted, a pre-processing step is performed in which a mask is applied to the extracted phase data to remove invalid noise information. According to Otsu's method, we used a global image thresholding [38]. This is required to increase the algorithm's processing speed and perform the post-processing steps described in section 3.2. In the following, $[\; ]$ represents the standard rounding function, and $\lfloor\;{} \rfloor $ symbolizes the floor function.

The 2π ambiguities are first determined using a simple linear extrapolation ${\gamma _{x,y}}(t )$ at time $t$:

In our context, the time difference $\Delta t = ({t - 1} )- ({t - 2} )$ between the two successive time points always maintains a value of one and can therefore be neglected. To determine whether a positive or negative 2π correction is required, the position of the extrapolation point ${\gamma _{x,y}}(t )$ and the wrapped phase ${\phi _{x,y}}(t )$ are evaluated:

Afterward, the integer ${q_{x,y}}(t )$ that minimizes the difference between ${\phi _{x,y}}(t )$ and the anticipated phase value ${\gamma _{x,y}}(t )$ is determined (see exemplary data in Fig. 3 under “Linear Extrapolation”):

Equation (1), ${q_{x,y}}(t )$ and ${j_{x,y}}(t )$ are used to calculate the unwrapped phase ${\varphi _{x,y}}(t )$.

In a subsequent assessment, the necessity for additional corrective actions is evaluated. A well-known strategy to investigate a signal is to analyze its first derivative, as it can reveal discontinuities or abrupt changes. Following analyses and corrective actions should be carried out if specific criteria are fulfilled:

If conditions 1, 2, and either 3 or 4 are fulfilled, supplementary steps are performed:

If the requirements are met, the sign of the first derivate $D_{x,y}^\varphi $ with $sgn:\,\mathrm{{\mathbb R}} \to \{{ - 1,\; 1} \}$ is initially evaluated to detect an abrupt direction change of the signal:

Figure 3 illustrates an exemplary data set under “Abrupt Direction Change Analysis” relevant to this step. The displacement curve displays an abrupt positive jump, which is also evident in the direction change of the first derivative. However, in the −2π scenario, this sudden change is no longer apparent, validating the suitability of applying a −2π correction in this specific case.

To enable comprehensive analysis, it is also advantageous to examine higher-order derivatives, as these provide additional critical information about the behavior of the signal and allow more complex movements to be examined. Therefore, the third-order partial derivative is also considered to identify less prominent phase ambiguities. A zero value implies that the derivate’s rate of change is constant, while the opposite is valid for a non-zero value. Considering the force of the air-jet together with our 2.45 kHz frame rate, corresponding to a time of 0.4 ms between two B-scans, we anticipate a steady force or acceleration, except during the air-jet's rising edge. As a result, a significant jerk in the phase signal is only likely during this period. However, the high frame rate should limit any observed jerk to a minimum. Assuming a mostly constant acceleration value, the displacement should experience minimal jerk over time:

Then, the impact of the estimated ${q_{x,y}}(t )$ value on the phase progression is examined using an environment-aware predictive approach. Since it is known that a wrong correction distorts the phase values of all the following ones, errors can be identified and rectified with this method. In the estimated case, a forecast is initially made for a phase value at the time $t + 1$. To utilize Eq. (1) again, it is necessary to calculate the parameters ${j_{x,y}}({t + 1} )$ and ${q_{x,y}}({t + 1} )$, involving a computation of ${\gamma _{x,y}}({t + 1} )$ using the steps outlined in Eqs. (4) to (6). These steps are then repeated for the time $t + 2$. Various tests have shown that two predicted phase values are sufficient for reliable error detection. However, adjusting the algorithm to different needs may involve considering varying numbers of forecasted phase values, requiring suitable adaptations. In addition, the effects of a ± 2π correction at time t are evaluated with this method. ${\xi ^{est}}$ indicates the predicted phase values for the estimated case. At the same time, ${\xi ^{ + 2\pi }}$ and ${\xi ^{ - 2\pi }}$ represents the phase response after performing a ± 2π correction on the estimated phase value, respectively (see example data in Fig. 3 under “Forecast Future Positions”). Next, it is analyzed whether the temporal evolution of the predicted phase signals $PEvo$ remains unchanged after performing a ± 2π correction. The results are compared with the initially determined phase direction:

Subsequently, the most probable phase progression case is determined. This is done by calculating the first derivative of the predicted phase values for all cases ${\dot{\xi }^{est/{+} 2\pi /{-} 2\pi }}$. Then, a median value is determined from all the pixels within ${\dot{\xi }^{est}}$. This median value is then compared with ${\dot{\xi }^{est/{+} 2\pi /{-} 2\pi }}$ at each pixel $({x,y} )$. This gives the absolute difference from the median $MD_{x,y}^{est/{+} 2\pi /{-} 2\pi }(t )$:

Furthermore, as previously described, an additional step is necessary to minimize the jerk in the phase signal (see Eq. (17)). Therefore, an error check based on the $MD$ values is conducted:

The signal is unwrapped by applying this process to each frame, and its true phase can be visualized over time (see Visualization 2, Visualization 3).

3.2 Post-processing

As the phase is highly sensitive and reliable results for the entire sample are needed for subsequent analysis and assessment, a post-processing step is required to remove any remaining random aberrations. For this, a local median filter detects and corrects single outliers. This involves determining the phase difference between vertically adjacent pixels in the last frame of the acquired OCT scan. If the deviation exceeds a threshold of π, it is corrected by a 5 × 5 sliding median window. The significant advantage of this approach is that only local outliers are corrected without losing structural information by smoothing the entire image. Visualization 2 illustrates that the displacement curve directly represents the force curve applied by the air-jet.

3.3 Stiffness calculations

Once the phase information is unwrapped, the sample deformation can be characterized, and the mechanical properties analyzed. One mechanical property often studied in this context is stiffness. It describes the degree to which a material deforms under a given load and is linked to its strength or elasticity. A value characterizing the stiffness ${k_{x,y}}(t )$ of the sample can be estimated by using the displacement $\Delta {d_{x,y}}(t )$ and the force ${F_{x,y}}(t )\; $ applied by the air-jet:

4. Robustness investigation of the algorithm

In this chapter, the robustness and functionality of the algorithm are comprehensively evaluated on two different tissue types. Except for section 4.4, all maps shown in this chapter represent the displacement of each pixel at the time of maximum load (or after approximately 22 ms; see exemplary force curve of the air-jet, Fig. 1) relative to the initial frame. Structures with low displacement are shown in red. In contrast, structures with high displacement are shown in blue. A rainbow color spectrum represents the transition from low to high, highlighting the wide range of variations.

4.1 Quantifying inhomogeneity detection performance

Various measurements were performed on an inhomogeneous PVA phantom to investigate whether different structures within the tissue can be detected. For this, the surrounding matrix with a Young's modulus of 820 Pa, the PVA inclusion with an elastic modulus of 2.26 kPa, and the transition between the two were measured with a peak force of 33 µN. Figure 4 shows the displacement maps, the 20 times averaged B-scans taken from the 3D volume, and the en face sample projections. In the center of the orange rectangle, the line scans were conducted.

 

Fig. 4. Evaluation of irregularity detection using inhomogeneous PVA phantoms. A composite image is shown consisting of the en face images, the 20 times averaged B-scans of the scan position taken from the 3D volume, and the corresponding displacement maps. Measurements were performed on the surrounding PVA matrix, the PVA inclusion, and the transition region between these two structures. The line scans were acquired in the center of the orange rectangle.

Download Full Size | PDF

The displacement maps show a striking contrast: the surrounding matrix has a larger displacement, represented by the blue color, while the PVA inclusion has a comparatively smaller indentation, represented by a yellow-greenish color in the displacement map. A similar color representation is observed when the transition from the surrounding matrix to the inclusion is measured. Therefore, the pseudo-colored displacement maps identify the different PVA structures. Since the results display the sample displacement obtained directly from the phase information, the force profile of the excitation source can be observed in the surrounding PVA matrix. However, the presented OCE imaging method and phase unwrapping algorithm enable effective visualization and identification of different mechanical properties of the tissue.

4.2 Reproducibility assessment of the measurements

After confirming the ability to visualize different structures within a tissue, the next step was to investigate the capability to obtain reproducible results. For this purpose, four measurements were performed in direct succession on an inhomogeneous PVA phantom and on an ex vivo WHO grade II meningioma brain tumor sample with an applied peak force between 30 and 32 µN. Figure 5 provides the applied force values specified for each measurement.

 

Fig. 5. Examination of the algorithm regarding the generation of reproducible results. Measurements were performed on a PVA phantom (first column) and on an ex vivo WHO grade II meningioma brain tumor sample (second and third column). The brain tumor sample was examined at two different positions. Four measurements were taken at each position in immediate succession. The corresponding displacement maps at the time of maximum load are shown. The applied force is indicated for each measurement.

Download Full Size | PDF

The first column illustrates the measurement results obtained with the PVA phantom. Following the methodology used in the previous chapter, the transition between the surrounding matrix and the PVA inclusion was measured and plotted. Referring to Fig. 5, a clear distinction is evident based on the color representation. Three additional measurements were taken at the same location to demonstrate the algorithm's reproducibility. In addition, four measurements were performed at two distinct positions of the ex vivo brain tumor sample. Again, a clear differentiation within the heterogeneous tissue can be observed. The difference between the first measurement and all subsequent ones was determined to quantify reproducibility. The mean difference and the standard error of the mean were then calculated. The determined variations were as follows: 249 ± 87 nm for the phantom, 132 ± 12 nm for meningioma brain tumor position 1, and 268 ± 7 nm for meningioma brain tumor position 2. Therefore, the repeated measurements show remarkably consistent results, indicating high reproducibility. It should be noted that slight variations may occur due to the successive measurements. Repeated measurements can influence the results, leading to minor deviations and changes. However, the results show that the algorithm can accurately determine and visually display displacement differences over multiple iterations.

Furthermore, the results were compared with the conventional MATLAB unwrap function to evaluate the algorithm's efficiency against a widely used standard method. The evaluation was performed on a WHO grade IV glioblastoma brain tumor sample with a peak force of 361 µN. The phase information was unwrapped using three approaches: the standard MATLAB function, the novel unwrapping algorithm without post-processing, and with post-processing. The results obtained in each case are shown in Fig. 6. The phase progression over the entire acquisition time was plotted for specific positions marked by black stars in the displacement map.

 

Fig. 6. Comparison between the advanced phase unwrapping and the standard MATLAB unwrap function. An elastography measurement of a WHO grade IV glioblastoma sample was conducted. The phase information was unwrapped using three methods: the common MATLAB unwrap function (A), the novel phase unwrapping algorithm without post-processing (B), and with post-processing (C). The phase progression was plotted for certain positions indicated by black stars in the map.

Download Full Size | PDF

The region of the strong slope is particularly noteworthy as most 2π discontinuities occur there. In this region, the conventional unwrapping function in MATLAB (A) has difficulties properly correcting the phase data. The poor performance of the standard MATLAB unwrap function is most evident in the fact that the displacement curve shows an upward motion, although the air-jet applies a downward force. In contrast, the novel unwrapping algorithm without post-processing (B) provides better results. However, certain parts of the phase map (Fig. 6(B) left) still exhibit errors, indicating that the correction was improper in these areas. This may be due to a low OCT signal in these pixels. After performing post-processing (C), these localized errors are corrected. Consequently, the presented unwrapping algorithm leads to a significant improvement in the accuracy of the phase correction. The advanced algorithm shows significantly better results, especially in areas characterized by fast and intense changes, where multiple 2π shifts occur within a short time.

4.3 Analysis of sample displacement with varying elastic moduli and forces

For validations, measurements were also conducted on homogeneous PVA phantoms with varying elastic moduli ranging from 875 Pa to 16.38 kPa and a mechanical load between 31 and 33 µN. The applied force values for each measurement are given in Fig. 7. The displacement maps of four exemplary phantoms are shown on the left side of Fig. 7. All maps are available in Supplement 1.

 

Fig. 7. Sample displacement study using homogenous PVA phantoms. Eight phantoms were evaluated, with the left showing the displacement maps of four representative phantoms with different Young's moduli. The displacement maps of all measured phantoms are available in Supplement 1. On the right side, a diagram shows the median displacement of the central A-scan at maximum load. The resulting sample displacement is shown in blue, and its reciprocal is in orange. The applied force values are listed for each measurement.

Download Full Size | PDF

It is evident from the figure that the determined displacement decreases with increasing elastic moduli. Furthermore, unlike the inhomogeneous datasets in subsection 4.1, the maps exhibit a nearly homogenous sample displacement over the entire scanning range. However, again, the force profile of the air-jet is noticeably visible, especially for the softer phantoms. The median displacement of the different phantoms was evaluated for the central A-scan and at maximum load to enable the data comparison without needing a force profile correction. The force $F = k \cdot x$ can be determined by the spring constant k and the displacement x. Assuming a constant force F, it can be derived that $k = 1/x$. The right side of Fig. 7 illustrates the results, with the acquired maximum displacement in blue and its reciprocal in orange. The blue curve shows an exponential decrease in displacement as the stiffness of the phantom increases. The observed linear trend of the reciprocal emphasizes that we operate within the linear range of the stress-strain curve, highlighting the significance of using PhS-OCT. Accordingly, apart from the requirement for a sufficiently high sample-dependent force, there is also the need for the lowest possible force, particularly in the case of highly viscoelastic and nonlinear brain tissue.

Moreover, the response of an ex vivo WHO grade I meningioma brain tumor to forces ranging from 4 to 243 µN was investigated, as presented in Fig. 8. In these measurements, tissue differences were challenging to detect below 40 µN. However, as shown in Fig. 8(A), higher structural variations become apparent at forces of 71 µN and above, resulting in greater contrast. At a specific pixel position, indicated by a white star in the maps, Fig. 8(B) illustrates the phase progression over the entire acquisition time. The graph indicates that the maximum indentation of the sample increases with a higher load, which corresponds very well with the theory of a material with linear stiffness. For a more detailed analysis, the displacement at the time of maximum indentation was plotted against the applied force, illustrated in Fig. 8(C) in gray.

 

Fig. 8. Investigation of repeated measurements on a WHO grade I meningioma brain tumor under varying excitation loads. A: The displacement maps illustrate the maximum indentation per pixel obtained at different forces between 4 and 243 µN. B: Phase progression over the entire acquisition time at a specific pixel position marked by a white star in the displacement maps in (A). Each color represents a different mechanical load. Artificial offsets were added for better visualization. C: Corresponding force-displacement curve at the time of maximum load (gray) and linear curve fit (red). The force applied is specified for every measurement.

Download Full Size | PDF

The results obtained reveal a noticeable non-linear trend curve. A linear fit (shown in red) was added to the graph for better clarification. It is worth noting that the curve shows a decreasing trend for forces above approximately 187 µN. The sample displacement curve should be distinguishable from the noise floor to determine mechanical parameters precisely. In Fig. 8(B), these features can be seen from a force of 71 µN. In Fig. 8(C), we remain in the linear regime up to an indentation of 3.61 µm for that specific sample type at the time of maximum load from the air-jet. Therefore, assuming a minimum signal-to-noise ratio of 10 dB, it should be ensured that the displacement noise remains at 300 nm or below. Since the required force depends strongly on the sample, this emphasizes the necessity of using appropriate, sufficiently high sample-dependent forces and a suitable phase correction algorithm. As reproducible results are obtained, the observed contrast can be reliably determined.

4.4 Reliability evaluation through histological correlation

For the final assessment, an ex vivo WHO grade I meningioma brain tumor was histologically examined to confirm the validity of the observed contrast at a mechanical stress of 87 µN. A photo and an en face image of the sample are shown in Fig. 9(A) and (B), respectively. The line scans were performed in the center of the orange rectangle, with the corresponding 20-times averaged intensity image taken from the 3D volume shown in Fig. 9(C). The respective H&E-stained histological image and stiffness map of the scan area are given in Fig. 9(D) and (E).

 

Fig. 9. Comparing histology and displacement map of a meningioma brain tumor. Photo (A) and en face image (B) of the sample. Line scans were obtained in the center of the orange rectangle. C: 20 times averaged OCT B-scan taken from the 3D volume. Corresponding histological section (D) and registered stiffness map (E). F - I: Enlarged views of (D) and (E) at the red and blue outlined rectangle.

Download Full Size | PDF

For initial tests, the stiffness was determined by dividing the sample displacement at the time of maximum indentation by the corresponding force of the air-jet, using Eq. (26). The sample displacement obtained directly from the phase information was used without applying any additional correction steps. Following normalization, an adaptive histogram equalization was performed on the stiffness data to enhance the contrast between the various structural components in the map. Thus, the stiffness information given here is regarded as pseudo values. Again, a rainbow color spectrum was used to improve the visual representation of the stiffness variations. Red represents high stiffness, and blue indicates low stiffness. As mentioned in subsection 2.1, a nonlinear image transformation was applied to the stiffness map to ensure precise overlay with the accompanying histological image (see Visualization 4).

In the H&E-stained image, dense and irregular connective tissue surrounds the tumor predominantly on its right side. They can be identified by their light pink staining, elongated structures, and low nucleus density. The connective tissue consists of collagen fibers and other proteins that support and stabilize the surrounding tissue. These characteristics give the tissue a higher stiffness and, as a result, greater resistance to deformation. The adjacent tissue is the tumor, distinguished by its darker staining and high nucleus density. It should be noted that the paraffin embedding of the tissue samples resulted in visible fissures and may have caused small sections to be missing from the surface.

The stiffness map identifies the connective tissue, which has a significantly higher stiffness than the tumor tissue in the upper left region. Furthermore, the map follows the diagonal orientation of the connective tissue. Figure 9(F)/(G) and (H)/(I) provide a magnified view of the two regions indicated by a red and a blue rectangle. The red-outlined region shows the boundary between the connective and tumor tissue. This boundary is also visible in the stiffness map, where the connective tissue appears stiffer in the lower, reddish region. Simultaneously, the tumor seems more elastic in the upper, bluish area. The enlarged region highlighted in blue marks another boundary between different tissue types. In this case, the entire region appears red and stiff since the tumor tissue is surrounded by connective tissue. Thus, the results show that the phase contrast generated by the algorithm can precisely and reliably detect structural variations within the tissue.

5. Summary and outlook

This study displays the performance and robustness of a novel one-dimensional unwrapping algorithm through extensive validation utilizing multiple PVA phantoms and ex vivo tumor samples. It allows accurate evaluation and analysis of sample displacements throughout the acquisition time, with high-resolution and pixel accuracy. This facilitates the calculation and evaluation of various mechanical properties. Depending on the excitation force and stiffness of the sample, accurate analysis of multiple 2π displacements is achievable. In the measurements presented in this paper, acquired using the MHz OCT imaging system, the maximum phase shift between two B-scans was always 4π, corresponding to a displacement of 0.96 µm. The presented results were verified by histological findings, demonstrating that the algorithm can achieve promising results. Therefore, the structural contrast from this research could be used for artificial intelligence (AI)-based brain tissue characterization algorithms.

However, although this remains largely unstudied, the algorithm was optimized for our MHz OCE setup. For effective jerk analysis, it is crucial to have a sufficient sampling rate to maintain a quasi-constant force. Slower sampling rates can impair the performance of the jerk analysis. Therefore, the algorithm may need adaptation for other, slower systems. Nevertheless, based on initial findings, our system could also unwrap the phase information by utilizing every third dataset, thereby reducing the processing time by a third. Moreover, to optimize processing speed, the detection interval can be limited to a specific segment within the 200 ms air pulse, for example, during the strong edge. However, comprehensive testing using multiple datasets would be required, and suitable algorithm adjustments should be implemented if needed. Depending on the system, the threshold value for the binary mask may also need to be adjusted.

The environment-aware predictive approach uses the median over the entire image for error detection. If smaller excitation areas are present, it may be necessary to use a more localized median. Furthermore, depending on the specific requirements, more than two predicted phase values can be considered during the unwrapping process to improve accuracy. However, such an extension would lead to longer processing times and necessitate additional algorithm changes. If any phase wraps remain after the first evaluation, the algorithm can be reapplied to the previously unwrapped data set. This iterative process can be used to detect and correct any remaining errors.

The primary objective is to adapt the system for use in the operating room. It aims to identify the boundaries between healthy and tumorous brain tissue and support the surgeon during tumor removal. Therefore, a suitable numerical motion correction strategy is desirable to apply the described method for in vivo imaging. Pre-registering phase datasets become valuable for intraoperative in vivo brain measurements. It allows for determining the unwrapped phase based on the registration outcomes, offering a potential solution to the problems caused by significant sample movement during unwrapping. Since entire volumes can be acquired quickly with MHz OCT, a complete 3D registration and motion correction is mainly straightforward. The objective is to mount the scan head on a robotic arm for utilization within the operating room. This integration will enable surgeons to scan various positions by controlling the arm's movements accordingly. However, potential angular forces from the air-jet on the sample need further investigation regarding their effects on phase correction and determining mechanical parameters. Using faster resonant scanners is also desirable to extend the scanning range beyond the current limit.

In addition, the force profile of the excitation source is visible in the outcomes. The Gaussian force profile, with its maximum in the center of the image, causes a non-uniform indentation of the sample over the entire scanned area. Consequently, lower displacements are often observed in the edge regions. Accordingly, it is crucial to determine and correct the force profile from all measured data. Some preliminary tests for further determination of the mechanical properties of the tissue and correction of the air-jet profile have already been demonstrated [39]. A more in-depth systematic evaluation will be published elsewhere.