CC BY
7
Optimization of Preliminary Vector Averaging for Improving Strain-Estimation Accuracy in Phase-Sensitive Optical Coherence Elastography
Alexey A. Zykov, Alexander L. Matveyev, Alexander A. Sovetsky, and Vladimir Y. Zaitsev*
A.V. Gaponov-Grekhov Institute of Applied Physics of the Russian Academy of Sciences, 46 Ulyanova str., Nizhny Novgorod 603950, Russia
*e-mail: [email protected]
Abstract. In this paper we present a method that significantly improves strain estimation quality in phase-sensitive optical coherence elastography (OCE), in which interframe strains are proportional to gradients of interframe phase variations. The analyzed phase-sensitive OCT signals can be conveniently represented as complex-valued vectors, the arguments of which correspond to the sought interframe phase-variations. To reduce the noise influence on the strain-estimation accuracy, these complex-valued vectors are usually averaged over a certain window. Increase in the window size can greatly reduce the noise. At the same time, for inappropriately big windows, this smoothing may cause strong distortion of the reconstructed spatial distribution of strain in comparison with the actual distribution, especially in regions of large gradients of interframe phase variations and, correspondingly, increased strain level. Therefore, automatic adaptive choice of the window size is proposed to maximize the noise reduction and at the same time to avoid possible distortion of spatial strain distribution in the reconstructed strain maps. © 2023 Journal of Biomedical Photonics & Engineering.
Keywords: optical coherence elastography; OCE; automatic choice of parameters; adaptive averaging.
Paper #9032 received 27 Oct 2023; revised manuscript received 24 Nov 2023; accepted for publication 24 Nov 2023; published online 26 Dec 2023. doi: 10.18287/JBPE23.09.040311.
1 Introduction
In recent years, besides straightforward perfection of obtaining structural images in optical coherence tomography (OCT), the further progress is strongly focused on the development of new OCT modalities/extensions [1]. Their realization requires deep analysis of series of OCT scans, for example, for realization of OCT-based angiography [2] and elastographic modality in OCT [3, 4]. In OCT-based elastography, starting from the first studies [5], special attention is paid not only to characterization of elasticity of biological tissues, but also in a broader sense the problems related to visualization of displacements and strains have been actively studied. In this context, especially in earlier works, as well as in some more recent studies (e.g., [6-9]), utilization of correlation-based methods of speckle tracking was considered as was
initially proposed in Ref. [5] by analogy with ultrasound elastography [10] and strain mapping in mechanical engineering [11]. In the recent decade significant attention is paid to the development of strain-estimation methods based on phase-resolved OCT [4, 12-15]. In contrast to correlation-based tracking, in phase-resolved approaches the scatterer displacements are conveniently visualized without the necessity of search operations by straightforward evaluation of interframe phase variations. Then spatially-resolved strains can be found by numerically estimating the phase-variation gradients.
For enhancing the quality of such OCT-based strain estimation, various improvements of finding interframe phase-variations and their gradients were proposed. For example, least square fitting of the phase-variation depth dependence combined with amplitude weighting was considered in Ref. [12] to reduce contributions of the most noisy small-amplitude pixels. In studies [16, 17],
optimization of time intervals between compared OCT scans was considered for estimation of slowly-varying interframe strains and finding the cumulative strains (for example, such slow-strain visualization is required for estimating deformations caused by insufficiently relaxed internal stresses when controlling the shape of prepared implants [18], for studying osmotic deformations [19] or polymerization processes [17]).
For estimation of "instantaneous" interframe strains, the vector approach proposed in Ref. [13, 14] has become widely used in recent years by various groups [20-22]. In this approach, phase-resolved OCT signals are represented as complex-valued vectors in the complex plane. The vector approach is very efficient computationally [4, 22] and opens convenient possibilities for various averaging and weighting procedures, parameters of which should be optimized taking into account particular features of the reconstructed strain distribution. Since the spatial structure of strains and resultant interframe phase variations may strongly vary during the acquisition of OCT-signals, as well as may be strongly inhomogeneous within the field of view, the optimal parameters of such averaging/weighting procedures should be correspondingly adapted.
In phase-sensitive OCE, interframe displacements of scatterers determine interframe phase variations, so that to find interframe strains one should estimate gradients of interframe phase variations. This estimation may be performed, e.g. by applying least-square procedures to the depth-dependence of interframe phase variations [12] or using the above-mentioned vector approach [13, 14]. For the vector method of strain estimation, recent study [23] considered optimization of such an important parameter as the choice of scale (or step) over which gradients of interframe phase variations are estimated. In the presence of measurement noises, the analyzed phasevariation distribution is characterized by random fluctuations, so that the choice of a too small step when estimating gradient (the minimal step is evidently equal to one pixel) may cause very strong errors in the gradient and strain estimation. Choosing a larger scale helps to reduce errors in the phase-gradient estimation, but too large gradient-estimation steps may also degrade the estimated strain if the chosen step exceeds the spatial interval, beyond which phase-wrapping occurs for the current level of local strain. Such multiple wrapping of interframe phase variations within an OCT scan is quite a typical situation for rather moderate strains ~10-3...10-2, for which the compared OCT scans are not yet strongly decorrelated and the interframe phase variations are well visible. In view of this, in work [23] a method was proposed for automatic choice of the differentiation scale corresponding to the local strain value. The main idea of the method is that this differentiation scale should be made so large that the estimated phase difference over the chosen scale approaches the maximal unambiguous value ±n rad (in other words, phase variation is close to wrapping over this spatial scale).
In the current paper we consider another important issue which also has to be optimized when finding displacement and strains using noisy phase-sensitive OCT signals. Namely, we will consider the problem of spatial averaging of interframe phase-variation maps to reduce noise level in the estimated strain. Generally speaking, the usefulness of averaging procedures for noise suppression in phase-sensitive OCT is fairly evident and in the vector method of strain estimation [14, 15] the averaging of complex-valued signals (also called phasors [24]) has been routinely used, but without adaptive optimization. However, similarly to the inappropriately large gradient-estimation scale, the useful signal may experience decrease in the resultant post-averaging amplitude because of phase-variation nonuniformity within a too large averaging window. On the contrary, for too small windows, the averaging does not yet appreciably help to reduce noises. Thus, similarly to the adaptive choice of gradient-estimation scale discussed in Ref. [23] for strain estimation using phasesensitive OCT signals, the size of averaging windows should also be adaptively chosen taking into account the spatial distribution of reconstructed strains.
It is important that in the problem of strain estimation, the averaging can be used at various stages of OCT-scan analysis. Namely, the averaging (which in what follows is called "preliminary averaging") can be applied to the initial results of pixel-to-pixel comparison of phasesensitive OCT data. When displacements of scatterers are spatially inhomogeneous (i.e., different from simple translational displacements) the corresponding interframe phase variations are also spatially inhomogeneous. Consequently, even if the sought strain is homogeneous, the interframe phase-variation is spatially inhomogeneous and for too large averaging windows the summed signals may become essentially out-of-phase. This effect is the main factor that limits the size of preliminary-averaging window.
Additional averaging procedures (which will be called "final averaging") can be applied at the following stages to the quantities containing phase-variation gradients proportional to the sought strains. These gradients usually exhibit much slower spatial variability and the window size for final averaging may be chosen significantly greater.
In the following sections we specifically consider the problem of adaptive automatic choice of the averaging-window size at the stage of preliminary averaging. Using both numerically simulated OCT scans and real OCT records we will demonstrate that such averaging may be rather efficient. Along with the adaptive choice of the phase-gradient-estimation scale, the adaptive choice of window size for preliminary averaging can be used to significantly improve the quality of reconstructed strain maps in phase-sensitive OCT. The result based on numerical simulations and processing of experimental OCT images will be complemented with analytical arguments elucidating the physical background of the considered averaging procedures.
2 Illustration How the Size of Preliminary-Averaging Window Affects the Quality of Interframe Phase Difference Estimation
Currently, the phase approach is most widely used in optical coherence elastography (OCE) when estimating displacements and strains [4]. An important advantage of the phase approach is that, unlike correlation methods, the axial displacement component field u can be found without the need for any search operations. The pixel amplitudes of phase-sensitive scans can be represented in the complex-valued form a(m, j) = A(m, j)exp[ip(m, j)] , where absolute
amplitude A(m, j) =| a(m, j) | and <p(m, j) is the phase. Axial displacements u (m, j) of scatterers lead to interframe variations in the phases 0(m, j) =p2(m, j)-p(m, j) , which can be directly found by pixel-to-pixel comparison of the two subsequent OCT images, such that
u(m, j) =
®(m,j) = arg[ a2 (m, j)a* (m,j)],
(1)
where the asterisk denotes complex conjugation.
Thus, displacement u(m, j) in every pixel with the
axial index m and lateral index j can be expressed via these phase variations as [12]:
¿0®(m, j) 4mi
(2)
where is the central wavelength of the OCT source and n is the refractive index.
It is assumed that the interframe displacements may be of supra-wavelength magnitude producing phase wrapping, but they should not be essentially supra-pixel (otherwise, the pixel-to-pixel comparison makes no sense). In real conditions, the acquired complex-valued
amplitudes al2(m,j) are characterized by some
measurement errors. Consequently, the complex-valued
product a2(m,j)a*(m,j) and phase-variation field
0(m, j) found using Eq. (2) are also noisy.
The amplitudes and phases of measurement noises (e.g., the reception-array noises of the OCT device) are independent for neighboring pixels, whereas the sought more regular interframe signal variations in contrast to measurement noises are slower varying over the OCT scan. Consequently, for suppression of noise level, the complex-valued scans a2 (m, j)a* (m, j) can be averaged over a window Wz by Wx pixels in size, where z and x correspond to the axial and lateral directions, respectively. Since quantity a2 (m, j)a* (m, j) represents a complex-valued vector in the complex plane, such averaging procedure can reasonably be called "vector averaging" [13, 14].
Fig. 1 Demonstration of influence of averaging window size on the resultant phase-variation maps using a numerically simulated example. (a) Interframe strain distribution adopted in the simulations, (b) numerically simulated interframe phase difference without preliminary averaging. Other panels show phase differences after preliminary averaging using various window sizes Wz x Wx in pixels: (c) 2*2, (d) 2*16, (e) 16*2, and (f) 16*16. Bigger windows, while being the best in one area, strongly distort the phase difference in others.
| o, 5 mm
a m
« m
I0
HSHPB- rr '
WB&MMm
mmmmm
Fig. 2 (a) Adopted for numerical simulations spatially inhomogeneous strain distribution characterized by three regions of axially homogeneous strains of 2T0-3, 4T0-3, and 8T0-3, respectively. (b) The corresponding interframe phase-variation map found without any averaging. The noise level in the simulated scans is SNR = 1.5 dB.
In real conditions the spatial variability of reconstructed strains is a priori unknown, so that in different directions and different scan parts the allowable sizes of the averaging window can strongly differ. Usually, to minimize the distorting influence of averaging, but still obtain some noise reduction, one has to use the minimal size of the averaging window 2x2 pixels, although in such a case the noise reduction may be insufficient, although in some scan regions significantly stronger noise reduction can be obtained.
To illustrate the importance of proper choice of the window size for preliminary averaging of complex-
valued quantity a2 (m, j)a* (m, j), in Fig. 1 we present a numerically simulated example assuming that the interframe phase variations are produced by the strain field shown in Fig 1a, where in various scan regions the strain is intentionally chosen either laterally or axially homogeneous. For the simulations shown in Fig. 1, OCT beam is assumed enabling depth-independent lateral resolution, details of the simulation approach are described in Refs. [25, 26]. Sufficiently strong, but quite realistic SNR=3.5 dB is adopted in simulated Fig. 1, the measurement noise being modeled by adding random complex-valued quantities to the complex-valued amplitude of every pixel, so that the directly found via Eq. (1) interframe phase variation is rather noisy as shown in Fig. 1(b). The other panels show the phasevariation maps obtained after averaging the quantity a2(m, j)a1(m, j) over windows with various sizes 2x2, 2x16, 16x2, and 16x16 pixels. Fig. 2(c) shows that the smallest averaging window 2x2 pixels does not yet distort the shapes of phase-variation isolines, but yields only moderate noise reduction. The other larger windows give much stronger noise reduction in some areas without distortion of isoline shapes, but in other parts of the scan the inappropriately chosen window sizes and shapes produce very strong distortions of the interframe phase variation. In these image parts, the phase-variation
gradients and, consequently, the reconstructed strains obviously will be strongly distorted.
These examples shown in Fig. 1 already qualitatively indicate that, similarly to the adaptive choice of gradient-estimation scale [23], the adaptive choice of the window size for preliminary vector-averaging is also important. To show the importance of the preliminary vector averaging for the strain-estimation quality in more direct quantitative way we consider another example with simpler strain distribution allowing for easier interpretation. Fig. 2(a) shows the assumed strain distribution in the form of three vertical stripes and Fig. 2(b) shows the corresponding interframe phase variation 0(m, j) defined by Eq. (2). The size of the scans is 160 pixels axially and 256 pixels laterally, with the axial and lateral distances between the pixel centers, hz=6 ^m and hx=10 ^m, respectively. As for Fig. 1, the illuminating beam has depth-independent radius w0=10 ^m with the central wavelength ^=1 ^m. The number of scatterers in the scan is 256000, which corresponds to typical density of cells in biotissues. The noise level in the simulated scans is rather high with SNR=1.5 dB.
Next, the three columns in Fig. 3 corresponding to the phase-variation map shown in Fig 2(b) demonstrate the reconstructed axial strains that are proportional to the axial phase-variation gradients estimated by the vector method [14] for various values of the axial inter-pixel distance g = 1, 2, and 4 pixels (in other words, various gradient-estimation scale the role of which is discussed in detail in Ref. [23]). The upper row in Fig. 3 shows the influence of the chosen gradient-estimation scale g without any additional utilization of preliminary vector averaging of quantity a2 (m, j)a* (m, j) directly from the phase variations shown in Fig. 2(a). The bottom row in Fig. 2 corresponds to the same gradient-estimation scales g = 1, 2, and 4 pixels, but with additional preliminary
vector averaging of quantity a2(m, j)a1*(m, j) over the window 4x4 pixels.
Fig. 3 Reconstructed strain maps based on phase-variation map shown in Fig. 2b for various values of gradient-estimation scale g and averaging-window sizes Wx = Wz of 1*1 (a-c) and 4*4 pixels (d-f). The first column is for g = 1, second for g = 2, and third column for g = 4.
Fig. 3 demonstrates that the quality of the reconstructed strain essentially depends on the processing parameters. It is also clear from this example that the choice of the gradient-estimation scale g earlier discussed in Ref. [23] and the window size Wx x Wz for
preliminary averaging of quantity a2(m,j)a**(m,j) before estimating the phase-variation gradient may give significant contributions to the improvement of the strain-reconstruction quality.
In the following sections we specifically consider how the window size for preliminary vector averaging affects the quality of strains reconstructed using the vector method [14]. This method is chosen because of such advantages as intrinsic high robustness with respect to various noises and high computational efficiency. The results will be presented in quantitative form based on both numerical simulations like for Fig. 2 and supported by analytical argumentation. Then the method of adaptive choice of the window size Wx x Wz for
preliminary averaging of quantity a2(m,j)a*(m,j) will be formulated, to complement the earlier obtained results [23] for adaptive choice of gradient-estimation scale g. To make the discussion of these problems clearer, firstly, the main steps of the vector method [23] will be recalled in the next section.
3 Brief Summary of the Vector Method of Strain Estimation
In the vector method of phase-gradient estimation pixelated OCT signals characterized by complex-valued amplitudes a(m, j) = A(m, j)exp[i -$(m, j)] till the last
stage are processed without explicit phase extraction, according to the main steps schematically shown in Fig. 4:
a) The first step is pixel-by-pixel multiplication of the deformed B-scans and complex conjugated reference one:
a2(m, j)a*(m, j) = b(m, j) = = B(m, j) exp[i • Ф(т, j)],
(3)
where Ф^, j) = cp2 (m, j) - p (m, j) and
B(m, j) = A2 (m, j) A1 (m, j) . The matrix b(m, j) contains interframe phase variations Ф^, j).
b) Second, to improve the signal-to-noise ratio (SNR), vector averaging of quantity b(m, j) is
performed using a small sliding window Wz x Wx in size:
b(m, j) = В (m, j) exp[/' • Ф (m, j)] =
Wz Wx
= X Xb(m'+ m, j'+ j),
(4)
m'=\ j'=\
The role of just this averaging step will be considered in detail further.
c) Then for finding the axial phase-variation gradients, the next step is finding the following matrix:
d(m, j) = D(m, j)exp[i •*¥ (m, j)] = = b (m + g, j)b *(m, j).
(5)
Fig. 4 The main steps of the vector method of strain estimation. All intermediate calculations are performed with complex-valued quantities, and the sought phase-variation gradient is extracted at the final stage.
It contains the phase difference ¥(m, j) between pixels with axial indices m and m + g, where the interpixel distance g is the scale, over which the phasevariation gradient is estimated. In the simplest case g = 1, however the use of a larger g > 1 (if there is no ambiguity in the interframe phase variations on this scale) can significantly improve the quality of strain estimation as was discussed in Ref. [23].
d) To further suppress noise in the reconstructed field of axial gradients of interframe phase variations, vector averaging of quantity d (m, j) can be used with a window of size M x J pixels:
d(m, j) = D(m, j) exp[j • T (m, j)] =
M J (6)
= Z Z d (m'+ m j'+ j)■
m'=1 j'=1
The size M x J of this window usually is significantly larger than in stage (b), because the scale, over which gradient ¥(m, j) noticeably changes, is usually much larger than the scale over which the interframe phase variation 0(m, j) noticeably changes.
e) Finally, the sought strain was expressed via the phase-variation gradient T (m, j) as:
Sz = T (m, j)/(2k0hzg),
(7)
where k0 = 2x/X0 is the central wavenumber of the OC-source spectrum.
In what follows step (b) will be called preliminary vector averaging and step (d) - final vector averaging.
4 Influence of Window Size for
Preliminary Vector Averaging on the Axial Strain Estimation: Numerical Evaluation and Analytical Arguments
To numerically evaluate how the vector averaging of
*
quantity b(m, j) = a2(m, j)a*(m, j) affects the quality of the strain estimation, we used a pair of simulated B-scans, for which the interframe strain adopted in the simulation and the corresponding phase difference &(m, j) are shown in Fig. 2. We recall that SNR for the simulated scans was equal to 1.5 dB and other simulation parameters were indicated in the discussion of Fig. 2.
Based on the simulated data, strain is calculated according to Eqs. (3)-(7). The strain is calculated for different averaging window sizes in Eq. (4). The gradient-estimation scale g in Eq. (5) is set to 4 and the final averaging in Eq. (6) for quantity d(m, j) (the argument T (m, j) of which is proportional to strain) is
performed using the final averaging window 16x16 pixels.
Then we calculate the dependence of the found mean strain value Smean and the standard deviation (STD) for strain on the preliminary averaging window size. In the discussed simulated example shown in Fig. 2 there are three regions with three levels of homogeneous strain, for which the isophase lines are laterally strictly horizontal, so that the results of averaging in the axial and lateral directions should essentially differ. To clearly show this difference, while incrementing the window size along one axis the window size along another axis was set to fixed value 2. The numerically obtained dependences are shown (in black lines) in Fig. 5 as a function of Wx and Fig. 6 as a function on Wz.
It is clear from Fig. 5 that for horizontally homogeneous interframe phase variation, the estimates of Smean in Fig. 5 remain close to the value adopted in simulations and the deviation STD(Wx) monotonically decreases with increasing Wx, which is intuitively expectable. In contrast, the dependences on the vertical scale Wz of preliminary vector averaging are rather nontrivial. For the estimated Smean( Wz) in Fig. 6 (upper row), a clear drop is visible for Wz > g (we recall that g = 4 for this example). The behavior of STD(Wz) (see Fig. 6, bottom row) is essentially non-monotonic with a minimal deviation for a certain Wz, for which its relation to g will be discussed in the following sections.
To elucidate the origin of these numerically revealed dependences we may present the following analytical argumentation. Let us first consider a case of positive homogeneous strain Sz that is independent of the depth, as well as of the lateral coordinate, so that the interframe phase difference increases linearly with depth as in the discussed examples. However, even in case of such horizontally uniform strain Sz we consider a somewhat more general situation, in which the displacements of scatterers Az may depend not only on depth z, but may be accompanied by a vertical translation shift dependent on x (such that the iso-phase lines are inclined rather than horizontal as in Fig. 2).
Fig. 5 Numerical evaluation of the influence of preliminary averaging window size Wx on the axial strain estimation for the horizontally homogeneous strains and the corresponding interframe phase variations shown in Fig. 2. The (a), (b), and (c) panels show the estimated strain values found for the fixed gradient-estimation scale g = 4, fixed vertical size Wz = 2 of the averaging window and variable horizontal size Wx. The (d), (e), and (f) panels show the standard deviation for the estimated strain as a function of Wx.
Fig. 6 Numerical evaluation of the influence of preliminary averaging window size Wz on the axial strain estimation for the horizontally homogeneous strains and the corresponding interframe phase variations shown in Fig. 2. The (a), (b), and (c) panels show the estimated strain values found for the fixed gradient-estimation scale g = 4, fixed horizontal size Wx = 2 of the averaging window and variable vertical size Wz. The (d), (e), and (f) panels show the standard deviation for the estimated strain as a function of Wz.
Thus, for such interframe displacements of the form Az = szz + sxx , the interframe phase difference is
Ф= 2k0Az = 2k0 (szz + sxx) . For pixilated scans this
phase difference can be written as
Ф(^, j) = 2k0 [(hzszm + hxsxj)] , where hz and hx are the
axial and lateral distances between neighboring pixel centers; quantities m and j are the axial and lateral pixel indices.
We also assume that some measurement noise is present and is represented by adding to each pixel a complex-valued random quantity characterized by a Gaussian distribution with zero mean and standard deviation a. For these assumptions, the registered pixilated OCT signal can be written as a(m, j) = S(m, j) + n(m, j) , where S(m, j) is the true
OCT signal without the noise and n(m, j) is the additive Gaussian noise characterized by a correlational function
(n(ml, jl)n*(m2, ji)) = °2$(m2 -m)s(j2 - j1) , where
8 is the Kronecker symbol. The complex-valued signal before deformation can be written as S1(m, j) = A ■ exp[ip(m, j)] , where the initial phases
<p(m, j) are determined by random positions of scatterers and the amplitudes are taken equal for simplification. After deformation it takes the form S2(m,j) = A■ expfa(m,j) + 1Ф(m,j)} . Performing all
the steps of strain calculation for the signal part when a(m, j) = S(m, j) using Eqs. (3)-(6) one obtains the
following averaged value d of quantity defined by Eq. (6):
d = A4sinc2 (k0szhzWz) ■ sinc2 (k0sxhxWx) ■
■ exp(i2k0 szhzg I
(8)
where the sinc-function in Eq. (8) is sinc (y) = sin (y) / y. As can be seen from Eq. (8), averaging over a window with sizes Wz and Wx leads to decrease in amplitude of | d | with zero value reached for k0szhzWz = n or
k0sxhxWx = n
Then the calculation of the mean strain estimate sz
and its standard deviation STD( sz) was performed based on Eqs. (3)-(7) in the presence of noise (i.e., for a(m, j) = S(m, j) + n(m, j)) with account for the above formulated noise properties and assuming sufficiently high output (i.e., after Eq. (6)) signal-to-noise ratio for estimating strain sz and its standard deviation STD( sz).
It appears that the mean value for estimated strain sz is unbiased and equal to the actual strain s if Wz < g . However, when Wz is greater than scale g , over which the phase-variation gradient is calculated, the performed
analysis indicated the appearance of noise-induced bias A(sz) in the estimation of strain s :
A(sz) = ¡0, for
Wz < g
-sz-j- 2 + —
A I A I sinc2
smc(2kQszhzg )
W
(9)
(k„szhzWz )sinc2 (k„shW ) Wz
for W > g.
These approximate analytical estimates are shown by red lines with empty circles in Figs. 5 and 6. The bias of the estimated strain for Wz > g described by Eq. (8) is clearly seen in the upper row of Fig. 6. It is clear that the analytical estimates well agree with the numerical results. The reason of the above-mentioned bias in the estimated
s is due to the fact that for too large averaging with Wz > g , the noise acquires some correlation within the gradient-estimation scale g. Therefore, Wz must be chosen smaller than the gradient-estimation parameter g, for which the adaptive choice is discussed in Ref. [23]. In view of the independence of strain of the lateral coordinate for the considered strain shown in Fig. 2, there is no similar limitation on the lateral size Wx of the window for preliminary averaging (see the lower row in Fig. 5).
For the same conditions and assumptions, calculations based on Eqs. (3)-(7) yield the following functional dependence for the standard deviation STD( s ) of strain in the presences of noise n(m, j):
STD(sz ) I 2 +-
(10)
Khzg4MJ ■ sinc1 (KszhzWz) ■ sinc2 (KsxhxWx)-WW
In this analytical estimate we focus on the functional dependences rather than exact estimation of pre-factor in this equation. For clearer comparison with numerical results this pre-factor can easily be adjusted. The corresponding superposition of numerical (black color)
and analytically found (red color) curves for STD( s ) is shown by red lines in the bottom rows in Figs. 5 and 6.
Again for lateral preliminary averaging the curves of the two types (bottom row in Fig. 5) show monotonic
decrease in STD( s ) approximately as l/y[Wx as follows from Eq.(10) for sx = 0 corresponding to Fig. 2 with strictly horizontal isophase lines. However, for inclined isophase lines (i.e. for sx ^ 0 ), function
sinc2(k0sxhxWx) in Eq. (10) indicates that, due to horizontally varying phase, the result of lateral vector averaging should become non-monotonic as a function of
Wx.
For the vertical averaging, the dependence of STD( s ) on the scale Wz of vertical preliminary-averaging is essentially non-monotonic (bottom row in Fig. 6) because of the vertically inhomogeneous phase
l
variation the axial gradient of which is determined by the axial strain level sz.
Eq. (10) makes it possible to predict the sizes Wz and Wx of the preliminary-averaging window, for which the standard deviation STD( S ) is minimal. It follows from Eq. (10) that minimal STD( S ) is attained in the maximum of function f (Wz ,Wx) = sine2 (k0szhzWz) • sine2 (k0SxhxWx) -JWWx that enters the denominator in Eq. (9). The minima of STD( s ) as function of Wz and Wx are attained independently of each other. Using the notation a = 2k0szhzW it is possible to find the maximum of
function f (a) = sinc2(a /2) •ja . The shape of this function shown in Fig. 7 indicates that its maximum occurs at axx/2 , therefore, STD( S ) has minimum
near this point. In what follows criterion a~n/2 will be used for determining the optimal window sizes for preliminary averaging in both directions.
Fig. 7 The shape of function f (a) = sine2 (a /2) •ja that determines the averaging-window sizes W and W , for which STD( Sz) reaches minimum. The latter corresponds to f (a) maximum that is reached around n/2.
Further, when choosing the optimal W one also should take into account the absence of bias for the estimated S , namely, according to Eq. (9) it should be W < g , where parameter g is the scale used in Eq. (5) to estimate the axial phase-variation gradient. Optimization of parameter g was considered in detail in paper [23]. Although the use of larger values g > 1 is favorable for enhancing the gradient-estimation accuracy, in the presence of phase wrapping (when numerous rainbow stripes are present in the strain-variation map within the visualized region as in Fig. 2(a) along the vertical direction), the axial phase difference for too big g experiences jumps of 2n rad. (phase wrapping) and, consequently, the estimated axial strain changes the sign. Using this criterion of the incorrect sign change for the
estimated strain, the condition of optimal gopl can be formulated as gopt <X0/(4szhz ) . Then, recalling that a = 2k0shzW and the optimal a^K/2 , we obtain that for W , the optimal value can be chosen as WT = gopt /2 = XJ(%szhz ).
If the isophase lines are inclined like in Fig. 1 in the presence of strain component * 0 , then the same condition a~K/2, similarly to W , can be used in the lateral direction and the optimal value for the horizontal window size can be obtained, W^' < X0 / (8sxhz ). If ^ = 0 and the interframe phase variation does not depend on the lateral coordinate (like shown in Fig. 2(a)), there is no such limitation on the optimal W , so that the only limiting condition on its choice is the desirable spatial resolution in the reconstructed strain maps.
5 Algorithm for Automatic Adaptive
Choice of Preliminary Averaging Window Size
The above presented analysis allows one to formulate a practical algorithm for the adaptive selection of window sizes WXzz for preliminary vector averaging given by Eq. (4).
The first step is calculation of vertical and horizontal phase gradients with Eqs. (3)-(6). Note that for lateral gradient T in Eq. (6) it must be calculated with a horizontal step g(x) instead of vertical step g(z). In terms of auxiliary parameter a considered in Fig. 7, the criterion of the optimal window size for preliminary averaging was formulated in the previous section as a = K /2. It was demonstrated in Section 4 that in the axial direction z , this criterion corresponds to condition Wzop' = g Op' / 2 , where g^p) is the optimal gradient-
estimation scale gop' in the axial direction. It was
discussed in detail in Ref. [23] that the optimal value of gop' for estimating phase-variation gradient corresponds
to such a scale g in Eq. (5), over which the phase
difference T(m, j )|g = tends to wrapping (i.e.,
experiences a jump by 2k rad). Notice, that in Eq.(6) phase difference T(m, j)\g is considered for the scale g
in the axial direction, but very similar arguments also relate to the lateral direction in case of laterally inhomogeneous interframe phase variations. Thus, for the lateral direction x , similarly to the axial one, the optimal window size should be found as WOp' = gOp / 2 .
After finding the so-defined optimal scales gOp) and gOp) for estimating phase-variation gradients in every image pixel (m, j), we obtain the mask of optimal sizes of preliminary-averaging windows in the axial and lateral directions, Wz x (m, j) = g(ptx) /2. If the so-found window
size is not integer, a closest smaller integer value must be taken.
The next step is the realization of preliminary vector averaging of quantity b(m, j) via Eq. (3) to find the
averaged matrix b(m, j) using the determined sizes
W^p (m, j) as described above. Next, quantity d(m, j)
containing the interframe phase-variation gradient ¥(m, j) can be estimated using Eq. (5) and the final averaging of d(m, j) is performed using Eq.(6) to obtain
d (m, j) and ¥(m, j).
Also using the so-found less noisy data b(m, j) one may iteratively repeat the above-mentioned procedures for estimating the scales g^f} [23] and, therefore, find
refined sizes of the averaging window Wop = g(xpi'] / 2,
which should yield improved estimates of the sought strain.
Examples of strain estimation with utilization of the above-described adaptive preliminary averaging are given in the following sections for both simulated and real OCT data.
6 Adaptive Preliminary Averaging: Simulated Examples
In the following examples we again used the approach described in details in Refs. [25, 26]. for simulating images for spectral OCT systems with "typical" parameters - a weakly focused illuminating beam (having a Gaussian profile with radius = 16 ^m weakly changing with depth), the central wavelength of 1.3 ^m (in the air) and a spectrum width of 90 ^m. Each A-scan consisting of 256 pixels covering 2 mm in depth (in air) is formed by superposition of 256 complex-values spectral components. The lateral spacing of A-scans is 16 ^m, and horizontally the B-scans consist of 200 A-scans. Discrete point-like scatterers were randomly distributed in the visualization area with an average distance between them of ~5 ^m, which corresponds to the typical cell density in biological tissues. The measurement noise was modeled by adding to each pixel random complex-valued number with Gaussian statistics. The noise level in the so-simulated scans corresponded to a signal-noise ratio (SNR) of 3.5 dB. For such fairly strong noise, the optimization of estimation of phase gradient and strain is particularly important.
Fig. 8 Color coded adaptively chosen window sizes Wx,z in the range 2...16 for preliminary averaging and illustration of the phase-variation map improvement. (a) Adaptively chosen Wz, (b) adaptively chosen Wx, (c) interframe phase difference for premilinary averaging with a fixed window 2*2 pixels, (d) interframe phase difference for preliminary averaging with adaptively chosen Wz and Wx shown in panels (a) and (b). The effect of noise reduction for adaptive averaging is clearly visible.
fixed Wx=2, Wz= 2, fixed Wx=2, Wz= 2, adaptive Wx, Wz g=3 adaptive g adaptiven
Fig. 9 Comparison of strains calculated for various combinations of averaging parameters and gradient-estimation scale g. The (a1) and (b1) panels are for fixed preliminary averaging window size Wz = Wx = 2 and fixed g = 3. The (a2) and (b2) panels are for fixed window Wz = Wx = 2 and adaptively chosen g. The (a3) and (b3) panels are for adaptive choice of all three parameters Wz, Wx, and g. (a1), (a2), and (a3) panels - unaveraged strains; (b1), (b2), and (b3) panels - averaged strains.
The simulated distribution of strains, according to which the interframe displacement of scatterers was performed, and the corresponding interframe phasevariation map are shown in Fig. 1. Such a strain distribution was chosen with an idea to present a large number of possible scenarios in one image to demonstrate the capabilities of the method. Strong horizontally homogeneous strains typical for a soft silicone calibration layer which is used for mapping the elastic properties of biological tissues [27, 28] were modeled in the top layer. The modeled weak homogeneous strains in the middle layer of Fig. 1a and surrounding intense rather complex horizontally nonuniform strains of positive and negative signs are also typical of some real experiments, e.g., related to visualization of thermomechanically induced strains in samples of eye cornea and cartilages [29] were also modeled.
There are several different characteristic areas in the phase difference in Fig. 1(b): dense in vertical direction and homogeneous in horizontal (top layer), slowly varying in vertical direction and homogeneous in horizontal (middle layer), almost vertical lines in the bottom middle part and inclined isophase lines in many other areas.
To the so-simulated data, the algorithm with adaptive preliminary averaging window size described in the previous section was applied. The resulting masks for chosen axial and lateral window sizes Wx and Wz are presented in Figs. 8(a) and 8(b), respectively. These
masks are not strictly regular. The irregularities are due to the fact that they are chosen for the noisy simulated data. Comparison of the use of the proposed optimized averaging method and averaging with a fixed window size of 2x2 pixels is shown in Figs. 8(c) and Fig. 8(d). The adaptive averaging shows significant reduction of noise.
Fig. 9 demonstrates comparison of the reconstructed strain maps obtained for several combinations of the gradient-estimation scale g and parameters Wx,z for preliminary averaging. The upper raw in Fig. 9 demonstrates the strain distributions found using only preliminary vector averaging of the interframe phase variation (given by Eq. (5)), whereas the lower row shows the strain maps after final vector averaging of quantity d (m, j) containing phase-variation gradients
(using Eq. (6)). The latter procedure is made for the same averaging-area size 12x12 pixels in all the cases. The left column in Fig. 9 is obtained with preliminary averaging using small fixed widow Wz = Wx = 2 pixels and fixed g = 3; the middle column is for the fixed window Wz = Wx = 2 and adaptively chosen g; and in the right column all three quantities Wz, Wx, and g are chosen adaptively.
It is important to emphasize that, since the window sizes Wz and Wx for adaptive preliminary averaging are closely linked to the adaptive choice of gradient-estimation scale g, in Fig. 9 we compare all 3 cases: (i) when neither Wz x Wx, nor g are optimized, (ii) when only g is optimized, and (iii) when both parameter g and
the size Wz x Wx of the preliminary averaging are optimized. The strain distribution obtained with all parameters chosen adaptively demonstrates appreciably better correspondence with the simulated strain distribution shown in Fig. 1(a). The main improvement is seen in the upper layer, where the strain is maximal the optimized gradient-estimation scale g cannot be made larger than 3 pixels and, therefore, improving in SNR can only be achieved by means of preliminary averaging over a larger window.
7 Adaptive Preliminary Averaging: Real Data Examples
This section describes the application of preliminary averaging with adaptively chosen parameters to process OCT-scans obtained in contact mode for mechanically-produced deformations of a three-layer structure consisting of two silicone layers and a cartilage sample placed between them similarly to studies described in Ref. [30]. The OCT system enabled B-scans of 2 mm in depth (in the air) and 4 mm in lateral direction. Examples of registered structural image and the interframe phase difference are presented in Fig. 10. Rapid axial variation of the interframe phase difference in the upper soft silicone layer in Fig. 10 is similar to the upper layer in the simulated phase difference in Fig. 1b. Cartilage is much stiffer than silicone, therefore, in cartilage the interframe phase difference varies slower with depth and can be averaged over larger areas.
In the following Fig. 11, results similar to the model Fig. 8 are presented. Figs. 11(a) and 11(b) are the adaptively chosen masks for axial and lateral averaging sizes Wz and Wx. Notice that the boundaries separating regions with different values Wz and Wx may look somewhat more irregular than the variations in the density of isolines in the phase map shown in Fig. 10(b). These irregularities are due to noises, the influence of which causes fluctuations in the estimated parameters
g(op and gOp (i.e., the estimated distances over which
the interframe phase-variation experiences wrapping in z- and x-directions). Correspondingly, the estimated
averaging sizes Wx,z = gOp'tZ / 2 also fluctuate.
Figs. 11(c) and 11(d) show the phase difference calculated with fixed averaging 2x2 and adaptive averaging. It is worth noting that although in the bottom part the masks (Figs. 11(a) and 11(b)) have irregular structure and might have lower values than optimal due to high level of noise, this fact does not spoil strain estimation compared to fixed 2x2 window but rather does not improve the results, which is not critical.
The comparison of strain maps similar to Fig. 9 for nonadaptive vs. adaptive averaging procedures is shown in Fig. 12. Improvement in the reduction of noise level is quite noticeable especially in comparison with strain shown in Fig. 12(a1) and found for fixed gradient-estimation scale g=3 and 2x2 window for preliminary averaging without final averaging via Eq. (6). While strain is completely masked by the noise in Fig. 12(a1), it becomes better visible for adaptive parameter g in Fig. 12(a2) and in Fig. 12(a3) where adaptive g and adaptive preliminary averaging are combined yet without final averaging via Eq. (6). The additional use of final averaging in the lower row of Fig. 12 shows noticeable quality improvement in comparison with the upper row. Thus, Fig. 12 clearly demonstrates that the adaptive choice of both the gradient-estimation scale g and the preliminary averaging window size is very important for improving the strain-mapping quality in OCT-elastography.
8 Discussion and Conclusions
This study for the first time considers an important question for strain reconstruction in phase-sensitive OCT elastography - adaptive preliminary vector averaging of interframe phase difference. Averaging of phase difference is usually performed over a small fixed window to resolve the smallest scale of interframe phase variations in phase-sensitive OCT scans.
B-scan Interframe phase difference
Fig. 10 OCT visualization of a sandwich structure consisting of a cartilage sample between two silicone layers. (a) Structural B-scan, (b) interframe phase difference corresponding to significantly larger deformations in the reference silicone than in cartilage.
Fig. 11 Panels (a) and (b) show adaptively chosen axial and lateral averaging-window sizes. Panels (c) and (d) show interframe phase difference calculated with fixed averaging window 2 *2 pixels (c) and adaptive averaging (d).
The use of a bigger preliminary-averaging window can greatly improve the strain estimation quality. However, such larger averaging windows can be used without reduction in the useful signal only in image regions with sufficiently smooth spatial variation of interframe phase difference. Therefore, adaptive automatic choice of the preliminary-window size is proposed. It is demonstrated the properly performed preliminary averaging may additionally noticeably enhance the quality of strain mapping.
The recent study [23] considered another important issue for mapping strains in phase-sensitive OCT - the adaptive choice of scale g, over which the interframe phase-variation gradients are estimated. In the present study we show that the adaptively chosen window sizes Wz and Wx for preliminary averaging are closely related to the optimal gradient-estimation scales, The maximal allowable size of
W =
" x, z
parameter g x, z is limited by the lateral and axial
variability of interframe phase variations, namely by the scales, over which wrapping of interframe variations occurs.
The results of Ref. [23] and the present study indicate that the quality of strain mapping initially increases <x g
and к ^WXz z . Notice that quite often the interframe
phase variations are significantly weaker dependent on the lateral coordinate than in the axial direction. Consequently, the averaging size Wx may be increased in comparison with Wz to additionally enhance strain mapping quality by preliminary averaging.
Adaptive selection of scales W and W for preliminary-averaging is performed for each direction separately using basically similar procedures to the choice of optimal gradient estimation scale in Ref. [23]. In both directions, the selection of Wx z is based on
estimating the scale, over which interframe phase variations experience wrapping (i.e., a jump by n rad.). We recall that this scale determines the optimal values of scales g0Pt for estimating gradients of interframe phase
variations [23]. In the previous sections we presented numerical and analytical argumentation that the optimal
averaging-window size Wx0,pzt should be chosen so that
change in the interframe phase variation is around n/2 in the lateral or axial directions, respectively, which means that W°pz' = g0p /2.
Here, it is important to emphasize a non-trivial result of the present study concerning the possibility to use only optimization of W with some fixed g or vice versa.
Fig. 12 Comparison of strains calculated using various averaging procedures. The (a1) and (b1) panels are for fixed preliminary averaging window size Wz=Wx=2 and fixed g=3. The (a2) and (b2) panels are for fixed window Wz=Wx=2 and adaptively chosen gradient-estimation scale g. The (a3) and (b3) panels are for the adaptive choice of all three parameters Wz, Wx, and g. The (a1), (a2), and (a3) panels are for the absence of final averaging via Eq. (6) and panels (b1), (b2), and (b3) panels with final averaging.
Namely, the gradient estimation scale gx,z may be chosen adaptively and used for some fixed preliminary-averaging window with a reasonably small size Wx,z (say, 2x2 pixels) as was considered in Ref. [23]. Certainly, in this case the noise reduction would not be as efficient as possible, but optimization of only gradient-estimation scales gx,z would not cause regular distortions in the reconstructed strains maps.
On the contrary, a not evident result of the present study is that it is not recommended to use a fixed gradient-estimation scale gx,z (say, 2 or 3 pixels) in combination with the adaptive choice of sizes W^p,' for preliminary averaging as described above. In this case, it is likely that the so-found averaging-window size WXOp,'
might significantly exceed the chosen fixed gradient-estimation scale gx,z. In this context, Eqs. (9) and (10) indicate that if the averaging-window size Wx,z exceeds scale gx,z, then the errors in the strain estimation may be even enhanced rather than reduced. Furthermore, the estimated strain acquires a regular bias that may be quite
significant (these effects are especially clearly seen in the right column in Fig. 6). However, these negative effects should not occur when adaptive optimization of the averaging-window sizes WXOpz' is performed
simultaneously with adaptive optimization of gradient-estimation scale gx,z.
Overall, the presented analysis and results of processing of numerically simulated and real OCT records illustrate the efficiency of optimization of the averaging procedures in combination with optimization of the gradient-estimation scale for improving strain-mapping quality in phase-sensitive OCT elastography.
Acknowledgments
The study was supported by the Russian Science Foundation grant No. 22-22-00952.
Disclosures
The authors declare that they have no conflict of interest.
References
1. R. A. Leitgeb, F. Placzek, E. A. Rank, L. Krainz, R. Haindl, Q. Li, M. Liu, M. Liu, A. Unterhuber, T. Schmoll, and W. Drexler, "Enhanced medical diagnosis for dOCTors: a perspective of optical coherence tomography," Journal of Biomedical Optics 26(10), 100601 (2021).
2. A. Zhang, Q. Zhang, C. L. Chen, and R. K. Wang, "Methods and algorithms for optical coherence tomography-based angiography: a review and comparison," Journal of Biomedical Optics 20(10), 100901 (2015).
3. K. V. Larin, D. D. Sampson, "Optical coherence elastography - OCT at work in tissue biomechanics," Biomedical Optics Express 8(2), 1172-1202 (2017).
4. V. Y. Zaitsev, A. L. Matveyev, L. A. Matveev, A. A. Sovetsky, M. S. Hepburn, A. Mowla, and B. F. Kennedy, "Strain and elasticity imaging in compression optical coherence elastography: The two-decade perspective and recent advances," Journal of Biophotonics 14(2), e202000257 (2021).
5. J. M. Schmitt, "OCT elastography: imaging microscopic deformation and strain of tissue," Optics Express 3(6), 199-211 (1998).
6. J. Rogowska, N. A. Patel, J. G. Fujimoto, and M. E. Brezinski, "Optical coherence tomographic elastography technique for measuring deformation and strain of atherosclerotic tissues," Heart 90(5), 556-562 (2004).
7. J. Rogowska, N. Patel, S. Plummer, and M. E. Brezinski, "Quantitative optical coherence tomographic elastography: method for assessing arterial mechanical properties," The British Journal of Radiology 79(945), 707-711 (2006).
8. A. Nahas, M. Bauer, S. Roux, and A. C. Boccara, "3D static elastography at the micrometer scale using Full Field OCT," Biomedical Optics Express 4(10), 2138-2149 (2013).
9. J. Fu, F. Pierron, and P. D. Ruiz, "Elastic stiffness characterization using three-dimensional full-field deformation obtained with optical coherence tomography and digital volume correlation," Journal of Biomedical Optics 18(12), 121512 (2013).
10. J. Ophir, I. Cespedes, H. Ponnekanti, Y. Yazdi, and X. Li, "Elastography: A quantitative method for imaging the elasticity of biological tissues," Ultrasonic Imaging 13(2), 111-134 (1991).
11. B. Pan, K. Qian, H. Xie, and A. Asundi, "Two-dimensional digital image correlation for in-plane displacement and strain measurement: a review," Measurement Science and Technology 20(6), 062001 (2009).
12. B. F. Kennedy, S. H. Koh, R. A. McLaughlin, K. M. Kennedy, P. R. T. Munro, and D. D. Sampson, "Strain estimation in phase-sensitive optical coherence elastography," Biomedical Optics Express 3(8), 1865-1879 (2012).
13. V. Y. Zaitsev, A. L. Matveyev, L. A. Matveev, G. V. Gelikonov, A. A. Sovetsky, and A. Vitkin, "Optimized phase gradient measurements and phase-amplitude interplay in optical coherence elastography," Journal of Biomedical Optics 21(11), 116005 (2016).
14. A. L. Matveyev, L. A. Matveev, A. A. Sovetsky, G. V. Gelikonov, A. A. Moiseev, and V. Y. Zaitsev, "Vector method for strain estimation in phase-sensitive optical coherence elastography," Laser Physics Letters 15(6), 065603 (2018).
15. V. Y. Zaitsev, S. Y. Ksenofontov, A. A. Sovetsky, A. L. Matveyev, L. A. Matveev, A. A. Zykov, and G. V. Gelikonov, "Real-Time Strain and Elasticity Imaging in Phase-Sensitive Optical Coherence Elastography Using a Computationally Efficient Realization of the Vector Method," Photonics 8(12), 527 (2021).
16. V. Y. Zaitsev, L. A. Matveev, A. L. Matveyev, A. A. Sovetsky, D. V. Shabanov, S. Y. Ksenofontov, G. V. Gelikonov, O. I. Baum, A. I. Omelchenko, A. V. Yuzhakov, and E. N. Sobol, "Optimization of phase-resolved optical coherence elastography for highly-sensitive monitoring of slow-rate strains," Laser Physics Letters 16(6), 065601 (2019).
17. Y. Bai, S. Cai, S. Xie, and B. Dong, "Adaptive incremental method for strain estimation in phase-sensitive optical coherence elastography," Optics Express 29(16), 25327 (2021).
18. Y. M. Alexandrovskaya, O. I. Baum, A. A. Sovetsky, A. L. Matveyev, L. A. Matveev, E. N. Sobol, and V. Y. Zaitsev, "Observation of internal stress relaxation in laser-reshaped cartilaginous implants using OCT-based strain mapping," Laser Physics Letters 17(8), 085603 (2020).
19. Y. Alexandrovskaya, O. Baum, A. Sovetsky, A. Matveyev, L. Matveev, E. Sobol, and V. Zaitsev, "Optical Coherence Elastography as a Tool for Studying Deformations in Biomaterials: Spatially-Resolved Osmotic Strain Dynamics in Cartilaginous Samples," Materials 15(3), 904 (2022).
20. S. Kling, H. Khodadadi, and O. Goksel, "Optical Coherence Elastography-Based Corneal Strain Imaging During Low-Amplitude Intraocular Pressure Modulation," Frontiers in Bioengineering and Biotechnology 7, 453 (2020).
21. M. Singh, A. Nair, S. R. Aglyamov, and K. V. Larin, "Compressional Optical Coherence Elastography of the Cornea," Photonics 8(4), 111 (2021).
22. J. Li, E. Pijewska, Q. Fang, M. Szkulmowski, and B. F. Kennedy, "Analysis of strain estimation methods in phasesensitive compression optical coherence elastography," Biomedical Optics Express 13(4), 2224 (2022).
23. A. A. Zykov, A. L. Matveyev, A. A. Sovetsky, L. A. Matveev, and V. Y. Zaitsev, "Vector method of strain estimation in OCT-elastography with adaptive choice of scale for estimating interframe phase-variation gradients," Laser Physics Letters 20(9), 095601 (2023).
24. B. Baumann, C. W. Merkle, R. A. Leitgeb, M. Augustin, A. Wartak, M. Pircher, and C. K. Hitzenberger, " Signal averaging improves signal-to-noise in OCT images: But which approach works best, and when," Biomedical Optics Express 10(11), 5755 (2019).
25. V. Y. Zaitsev, L. A. Matveev, A. L. Matveyev, G. V. Gelikonov, and V. M. Gelikonov, " A model for simulating speckle-pattern evolution based on close to reality procedures used in spectral-domain OCT," Laser Physics Letters 11(10), 105601 (2014).
26. A. Zykov, A. Matveyev, L. Matveev, A. Sovetsky, and V. Zaitsev, "Flexible Computationally Efficient Platform for Simulating Scan Formation in Optical Coherence Tomography with Accounting for Arbitrary Motions of Scatterers," Journal of Biomedical Photonics & Engineering 7(1), 010304 (2021).
27. A. A. Sovetsky, A. L. Matveyev, L. A. Matveev, E. V. Gubarkova, A. A. Plekhanov, M. A. Sirotkina, N. D. Gladkova, and V. Y. Zaitsev, "Full-optical method of local stress standardization to exclude nonlinearity-related ambiguity of elasticity estimation in compressional optical coherence elastography," Laser Physics Letters 17(6), 065601 (2020).
28. M. A. Sirotkina, E. V. Gubarkova, A. A. Plekhanov, A. A. Sovetsky, V. V. Elagin, A. L. Matveyev, L. A. Matveev, S. S. Kuznetsov, E. V. Zagaynova, N. D. Gladkova, and V. Y. Zaitsev, "In vivo assessment of functional and morphological alterations in tumors under treatment using OCT-angiography combined with OCT-elastography," Biomedical Optics Express 11(3), 1365-1382 (2020).
29. O. I. Baum, V. Y. Zaitsev, A. V. Yuzhakov, A. P. Sviridov, M. L. Novikova, A. L. Matveyev, L. A. Matveev, A. A. Sovetsky, and E. N. Sobol, "Interplay of temperature, thermal-stresses and strains in laser-assisted modification of collagenous tissues: Speckle-contrast and OCT-based studies," Journal of Biophotonics 13(1), e20190019 (2020).
30. V. Y. Zaitsev, A. L. Matveyev, L. A. Matveev, G. V. Gelikonov, O. I. Baum, A. I. Omelchenko, D. V. Shabanov, A. A. Sovetsky, A. V. Yuzhakov, A. A. Fedorov, V. I. Siplivy, A. V. Bolshunov, and E. N. Sobol, "Revealing structural modifications in thermomechanical reshaping of collagenous tissues using optical coherence elastography," Journal of Biophotonics 12(3), e201800250 (2019).
