CC BY
4
Numerical Study of Supra-Wavelength Axial Motion Compensation in Contact-Mode Optical Coherence Angiography Using Fourier-Shift Procedures
Alexey A. Zykov*, Alexander L. Matveyev, Lev A. Matveev, and Vladimir Y. Zaitsev
Institute of Applied Physics Russian Academy of Sciences, 46 Ulyanova str., Nizhny Novgorod 603950, Russia
* e-mail: [email protected]
Abstract. Label-free angiographic methods based on optical coherence tomography (OCT) visualize blood vessels utilizing detection of red blood cells motion against surrounding static tissue. However, in practice, the surrounding tissue is never still due to natural motions of living organisms (e.g., breathing or heart beating). To mitigate large scale motions of the tissue relatively to the OCT probe the tissue examination can be performed in contact mode. In such a case, however, the OCT probe inevitably exerts some pressure onto the tissue, so that bulk motions lead to interframe deformations and depth-dependent tissue displacements, which have to be numerically compensated prior to angiographic visualization. Usually, sufficiently small deformations primarily affect pixel phases in OCT images rather than pixel amplitudes, and, therefore, phase-only compensation of the masking motions may be fairly sufficient. However, in case of larger strains and supra-wavelength displacements, larger inter-scan phase variations of the order of several periods lead to the appearance of pronounced "decorrelation noise" in which variations in pixel amplitudes and phases are combined. This effect significantly degrades the quality of the final OCT-angiography images. In this paper, we present a new method allowing to a significant degree to compensate this phase-amplitude decorrelation caused by spatially-inhomogeneous supra-wavelengths displacements. This compensation is based on the Fourier-shift theorem, which allows one to back-shift fragments of the deformed OCT-scans to their initial positions before deformation. At the same time variations of pixels due to the motion of blood particles within smaller-in-size vessel cross sections are retained. Although such backshifts do not compensate relative motions of sub-resolution particles, this procedure efficiently reduces decorrelation even for fairly big spatially-inhomogeneous displacements and leads to much lower signal variability outside blood vessels while preserving high variability inside. The proposed compensation method is compared to the earlier proposed phase-only compensation using simulated data. Pronouncedly lower strain-induced artefacts and much higher contrast between blood vessels and background are demonstrated. © 2022 Journal of Biomedical Photonics & Engineering.
Keywords: optical coherence angiography; OCTA; bulk motion compensation; microvasculature visualization; blood vessels.
Paper #3551 received 13 Oct 2022; revised manuscript received 11 Nov 2022; accepted for publication 15 Nov 2022; published online 8 Dec 2022. doi: 10.18287/JBPE22.08.040303.
1 Introduction
Optical coherence tomography (OCT) is one of the most innovative and rapidly evolving techniques in recent years. Beyond improvement of structural visualization, including such problems as increasing imaging speed and resolution, other OCT-based modalities are evolving. One can mention polarization-sensitive OCT imaging [1, 2], speckle contrast visualization [3], OCT elastography - mapping of Young's modulus of a tissue [4-7], and the most broadly used OCT angiography (OCTA) - visualization of blood vessels [8-11].
OCT-based angiography enables label-free imaging of blood microcirculation based on singling out motion of red blood cells (RBCs) against surrounding "solid" motionless tissue. This is an important distinction of OCTA from such angiographic methods as opto-acoustic imaging based on selective absorption of optical waves with certain wavelengths [12]. The latter effect is independent of RBCs motion, so that perfused vessels with moving blood and the ones with blood flow stasis are not easily differentiable.
For OCTA, differentiation of blood flow from surrounding "motionless" tissue is based on detection of temporal speckle variability in OCT images in regions of blood vessels. Erythrocytes are the main scatterers inside blood vessels and it is their motion which leads to increased temporal speckle variability in vessel cross sections. Particular OCTA techniques use various approaches to the variability estimation. Historically, first were developed methods utilizing the Doppler frequency shift [13-15]. Some methods are using only variations in amplitude/intensity (speckle variance) [8] or only variations of phase (phase variance) [16] for a series of B-scans obtained for the same position. One should also mention approaches as differential operations [17] and evaluation of temporal decorrelation for a series of sequential OCT scans acquired for the same position [18, 19].
An important point is that often OCTA methods use acquisition of a pack of B-scans from the same position with a subsequent step-wise shift to the next position. Suchjump-like steps may lead to transitional oscillations. To eliminate this effect, a scanning pattern with a smooth shift from B-scan to B-scan with 10-14 fold self-overlapping in the direction of slow scanning was used in this study similarly to work [20]. In this approach "solid" biological tissue also exhibits some variability from B-scan to B-scan but this variability is much lower than inside blood vessels, therefore the latter ones still can be singled out.
The OCTA method used in this work is high-pass filtering (detailed discussion can be found in Ref. [20]), which naturally utilizes the described above smooth scanning protocol. Briefly, in this method, B-scans are high-pass filtered in the slow scanning direction over which the B-scans are significantly self-overlapped. "Static" tissue does not change much from B-scan to B-scan. This fairly slow variability corresponds to low frequencies in the spectrum whereas blood vessels have
high variability, and therefore correspond to high frequencies. By suppressing low frequencies (which basically corresponds to high-pass filtering) one suppresses the static tissue. As only a fairly small portion (~10) of B-scans are overlapped and it does not make sense to compare B-scans that do not overlap, high-pass filtering can be performed using a sliding window. Furthermore, the Fourier transform operations for finding the spectrum, high-pass filtering and the subsequent inverse Fourier transform for obtaining again the timedomain signal can be realized via convolution of the acquired sequence of time-domain signals with a properly chosen temporal sliding window. Its length is of the order of the number of overlapping B-scans (in the below-presented examples this number is chosen L = 7) and the Fourier spectrum of the sliding convolution window corresponds to the form of the spectral window required for high-pass filtering [20].
Actually for all OCTA methods based on detection of increased temporal variability of signal from the regions of vessels, a very serious problem is that the surrounding "static" living tissue is not static at all and usually exhibits fairly intense physiological motions (due to breathing, heart beating, etc.) These motions often may strongly mask the sought signal variability related to the blood circulation, which significantly corrupts the quality of OCTA images if no additional measures are undertaken to either physically immobilize the "solid" tissue or to digitally compensate them. Physical tissue immobilization is often used in model studies on animals (in particular, special immobilizing dorsal window chambers can be used like, e.g., in study [21]), however, immobilizing devices often are not suitable even in experiments with animals and, moreover, cannot be used for examining patients. Alternatively, various methods of numerical compensation of masking tissue motions have been considered (see, e.g., review [22]).
In particular, for non-contact OCT examinations, the physiological motions of "solid" tissue cause predominantly translational shifts of compared OCT scans. In methods based on comparison of complex-valued OCT signals (i.e., the signals with amplitude and phase) the signal variability is dominated by axial motions rather than lateral ones (similarly to the domination of axial Doppler effect over the transversal Doppler effect). Axial motions with essentially sub-wavelength amplitudes mostly affect the signal phase and can be efficiently compensated using estimation of depth-averaged phases of compared in-depth A-scans [23]. However, motions with significantly stronger (supra-wavelength) axial displacements cause significant phase variations of the order of several periods and, moreover, are accompanied by appreciable variations in the signal amplitude. In combination they lead to high "decorrelation noise" which strongly masks the sought variability caused by blood-particle motions and, correspondingly, causes strong degradation of the final OCT-angiography images.
To prevent such large-amplitude motions typical of non-contact examination of living tissues, for situations
in which the use of tissue-immobilizing devices is impossible (especially when working with patients), OCTA can be realized in contact mode as proposed in Ref. [20]. However, along with suppression of large-amplitude translational motions in contact mode, the pressure of the OCT probe onto the tissue causes deformations (strains) in the tissue. The resultant strain-induced displacements of scatterers still may be fairly intense (their amplitude may be essentially supra-wavelength). Next, for strain-induced deformations in contact mode, displacements of scatterers are essentially spatially inhomogeneous unlike translational displacements in non-contact mode. Despite smaller amplitude of displacements in comparison with "typical" displacements in non-contact mode the above-mentioned strain-related interframe deformations still produce a rather strong masking effect and may significantly reduce the quality of angiographic images. Thus, such strain-induced displacements also need to be compensated. However, as demonstrated in Ref. [20], in contrast to translational motions, this compensation cannot be depth-averaged along the entire scan depth. Rather, the compensation procedure should take into account the spatial inhomogeneity of the displacements. Paper [20] demonstrated a variant of such compensation which was based on signal-phase correction (similarly to the translational case [23], but with additional accounting for depth-dependence). Compensation in such a form, nevertheless, retains to an appreciable degree the deformation-induced masking decorrelation of the scans in areas of "solid" tissue, especially in regions of supra pixel displacements, when the signal phase and amplitude both contribute to the scan decorrelation.
In this paper, we consider an improved compensation method to mitigate masking effects in the case of such spatially-inhomogeneous deformation-induced masking motions of the "solid" tissue. The method utilizes a spectral approach to fairly accurately locally back-shift fragments of deformed complex-valued OCT scans (including supra-wavelength displacements). To demonstrate the efficiency of the proposed approach and compare it with the earlier used pure phase correction, we utilize processing of 3D data sets of simulated OCT data using the flexible simulation platform reported in the earlier paper [24]. The utilization of simulated data sets is very convenient because unlike physical experiments parameters of both the masking tissue motion and the sought motion of scatterers mimicking erythrocytes are precisely controllable and can be very flexibly varied. Furthermore, measurement noises with controllable level can also be easily simulated. For such highly controllable conditions enabled by the used simulation we demonstrate that the proposed approach enables significant improvement of angiographic OCT-based visualization in contact mode.
2 Simulation Model
There are different approaches to simulation of OCT scans. The most known are Monte-Carlo methods [25-27] in which pseudo-random generation of
trajectories for a huge amount of photons is made to calculate their propagation and scattering. It performs really well in a situation when a single scan or a few scans are required. However, due to high computational demands, the Monte-Carlo approach is not very suitable for generation of large blocks of B-scans with moving particles to imitate OCT-based visualization of blood flow. In this study, we use the description of OCT-scan formation by weakly focused illuminating beams presented in Ref. [24], in which elements of the earlier proposed simplified model [28] (for one-dimensional illuminating beams with uniform cross section) and model [29] with rigorous accounting of illuminating-beam focusing are combined. The model form [24] is computationally rather efficient and at the same time fairly realistic for describing scans for widely used OCT systems with weakly focused beams that enable approximately depth-independent lateral resolution. In this approximation, it is possible to account for a Gaussian distribution of the beam amplitude in the lateral direction, whereas the phase front curvature can be neglected, which was verified by comparison with more rigorous model [29]. Next, the used model utilizes the commonly used approximation of ballistic scattering from a set of localized discrete scatterers [30]. Initially these point-like scatterers are seeded within the 3D imaged volume. Their positions can be chosen arbitrary, so that for imitating real tissues it is reasonable to randomly distribute such discrete scatterers with the density corresponding to the typical concentration of biological cells (approximately 5 ^m between neighboring randomly located scatterers). The scattering strength for the scatterers can be also varied (for example, to imitate regions with various brightness). Certainly, scattering strength with random distribution around some mean value can easily be realized, but it was verified that even for identical scatterers, the speckle statistics becomes "fully developed" (with speckle amplitudes fairly close to Rayleigh distribution) starting from a concentration higher than 1.5-2 particles per sample volume [28]. The Gaussian profile of the scanning illuminating beam readily allows one to describe how the scatterers gradually leave and enter the beam cross section. The scatterers can experience arbitrary displacements between the consequently formed OCT scans. Due to this feature, the model can easily simulate both the blood flow and various masking motions of the surrounding tissue [24].
Other main features of spectral-domain OCT scanners (the central wavelength, the number of received discrete spectral components and overall shape of the illuminating-beam spectrum) can also be readily accounted for the used model [24]. These spectral characteristics determine the resultant axial resolution (axial size of the point-spread function in the OCT image), the total imaging depths limited by the step between neighboring spectral components and the number of pixels in the formed complex-valued A-scans. Another feature that is important for the present study is that the used simulation method readily allows one to
imitate various measurement noises with any required signal-to-noise ratio by adding to every pixel random complex-valued quantities.
With all these features, using the 4-core (Intel Core i5) laptop Lenovo IdeaPad S340 the model generates 2D image 250 by 40 pixels with 8000 point-like scatters in 7 sec and 3D image 250 by 25 by 90 pixels with 170000 scatterers in a reasonable amount of time - 4 h. For the chosen scan sizes, this amount of scatterers imitates the typical density of living cells. It is worth to mention that simulation time can readily be reduced by using more computationally powerful facilities. The other parameters for simulation of the OCT images were chosen corresponding to the parameters of a typical experimental setup: central wavelength of the spectrum -1.3 mm in the air, spectrum width ~100 nm and the radius of the illuminating beam is ~15 ^m. The pixel size is 4 ^m in the axial z-direction and 16 ^m in the fast scanning x-direction in the B-scan plane. In the slow-scanning y-direction the inter-B-scan step is 2.14 ^m, this step being ~14 times smaller than the beam diameter.
3 Bulk Motion Compensation
As was mentioned above, a great challenge of angiographic visualization is that tissues surrounding blood vessels are not static due to natural motions of living objects such as heartbeat, breathing, etc. Since in OCTA the vessels visualization is based on estimation of signal variability, one must mitigate the masking effect of tissue motions when performing angiographic signal processing. There are two main options to do this.
The first option is to physically immobilize the examined tissue, e.g., using a dorsal chamber with a transparent window which is often used in experiments with animals [21, 23]. Such immobilization helps to deal with large-amplitude translational tissue motions relatively to the OCT probe. However, small-amplitude residual motions may still occur and their masking effect is mostly related to phase variations of the OCT-signal due to the axial component of these small-amplitude motions. Indeed, phase of the pixels in OCT scans is much more sensitive to axial motions rather than lateral, and small sub-wavelength masking motions mostly lead to phase variations. For the discussed translational motions, the resultant inter-scan phase variation is the same across all the depth. These fairly small depth-independent phase variations can be compensated numerically by calculating and compensating phase
difference Ap of the compared in-depth A-scans with complex-valued amplitudes A\ and A2; the phase difference is averaged along the entire depth:
Ap = aig I E A (z) 4 (z)
(1)
where z is a pixel index, N is the number of pixels in an axial direction, Ai is the reference A-scan and A 2 is the compared A-scan. Thus, the corrected complex-valued
amplitude A2 is acquired by the correcting phase factor
as follows:
4 = Â2 exp(iAcp).
(2)
This procedure helps to efficiently match the interscan phases in the regions of tissue around the vessels. On the contrary, in the cross-sections of vessels, the motions of scatterers are essentially independent of the bulk tissue motions, so that such depth-averaged inter-A-scan phase correction still preserves the phase variability within vessel cross sections. Then, this local signal variability can be singled out in one or another way (e.g., by estimating the signal variance [21] or applying highpass filtering [23]). The idea of the compensation method is schematically shown in Fig. 1.
The above-discussed immobilizing devices usually are not acceptable for application on patients. For this reason, an alternative approach was proposed in Ref. [20] to perform angiographic visualization in contact mode (when an OCT probe is in direct contact with the examined tissue, which may be quite acceptable for both animals and patients). The direct contact allows one to prevent large-scale tissue motions relatively to the OCT probe, however, the pressure of the OCT probe exerted on the tissue may cause tissue deformation instead of translational shift in non-contact mode. The deformation-produced displacements of scatterers in the tissue are not spatially homogeneous, and therefore phase compensation should be used in a depth-dependent manner using a sliding averaging window [20]:
A p( z, x) = arg
( z +Zp-1 x+Xp-1 ^
E E Bi(z' x')
z '=z x - x
V'B*(z' x')
(3)
A
Fig. 1 Schematically shown bulk motion compensation in case of non-contact mode. Red rectangles represent A-scans and blue circles - blood vessel. Depth-averaged phase is adjusted for subsequent A-scans, but for the local motions within the vessel cross-section the phase variations significantly differ from the averaged phase-variation values, such that they are preserved and can be singled out using any angiographic method.
Initial structural B-scan B-scan after deformation Initial phase difference Compensated phase difference
Fig. 2 Simulated OCT images: (a) structural B-scan before deformation, (b) B-scan after deformation (axial strain is 7.5*10-4). Both B-scans are very similar, so that blood vessels cross sections are not visually distinguishable in the structural images. (c) initial pixel-to-pixel inter-scan phase difference; (d) phase difference given by Eq. (5) after the described phase compensation. In both pictures three blood vessel cross sections appear as noisy regions due to the random motions of scatterers in the blood flow. The scale bar is 0.2 mm in all panels.
Such that the corrected complex-valued amplitude of the compared B-scan takes the form:
л
B2 (z, x) = B2 (z, x) exp(/'Д ф( z, x)), (4)
where z and x are indexes of a pixel in axial and lateral directions, respectively; Zp and Xp are sizes of the sliding window in pixels; B\ is the reference B-scan, B2 is the deformed B-scan, B2 is B2 after phase correction. The hat symbol л over phase variation Дф indicates that it is averaged over the sliding window. The averaging window size should be larger than the characteristic scale of the point-spread function, but at least several times smaller than the characteristic scale of inhomogeneity of the deformation-induced displacements. For visualization purpose we calculate the phase difference between the reference and compensated B-scans:
Дф'( z, x) = arg (B2 (z, x) B* (z, x)). (5)
Results obtained with this compensation method are presented in Figs. 2c and 2d showing the phase difference between a pair of simulated B-scans passing through cross sections of three vessels, in which the moving particles experience inter-scan displacements to imitate blood flow.
An important point about the compensation of bulk motion should be mentioned. For translational masking motions (as in Ref. [23]), the phase variations for all compared scans can be compensated towards the first reference scan. However, for spatially-inhomogeneous strain-induced displacements, strain may monotonically accumulate in time, so that the phase difference between sufficiently distant scans may become rather large and even experience wrapping on a scale of a single pixel, so that the phase-difference distribution may become rather noisy and cannot be properly compensated. However, it should be recalled that for gradually shifted position of B-scans planes, the phase comparison actually makes sense only for a limited portion of B-scans that are
partially overlapped. For even larger separation, the B-scans visualize quite different tissue portions and bear independent information. Therefore, for singling out vessel cross sections with moving scatterers, the informative data are given by comparison of fairly small portions of partially overlapped scans. Therefore, it makes sense to compare phases only within such limited scan portions and perform inter-scan phase compensation to the central B-scan of every self-overlapping small pack. Both variants of choosing scan portions for phasevariation compensations are schematically shown in Fig 3.
n+6
^ self-overlap
(a) (b)
Fig. 3 Illustration of different approaches to phase compensation in terms of chosen reference B-scan (colored in yellow). In panel (a) all the B-scans are phase-compensated toward the first B-scan, in this case, OCTA image may suffer from cumulative strain; in (b) B-scans are compensated only in a small pack which is used for convolution with filtering function. This sliding phasevariation compensation within self-overlapping scan packs eliminates the strain-cumulation issue, but the is more computationally demanding.
Each of such packs of self-overlapping scans can be subjected to high-pass filtering, which can be conveniently made by performing convolution with a filtering function sliding along the entire pack as proposed in Ref. [20]. In this case, eventual cumulation of strain and possible phase-estimation errors are much smaller, which improves the quality of phase compensation and final OCTA image.
This improvement for such sliding phase compensation is attained at the expense of a larger number of phase-comparison operations because roughly
N times more operations are required, where N is the length of a filtering window that covers portions of self-overlapping scans. Now, after the discussion of the different approaches to compensation of the bulk-motion induced phase-variability along the time axis (illustrated in Fig. 3), we return again to the features of interframe compensation of depth-dependent motions illustrated in Fig. 2. This method is based on only phase compensation with a window sliding along the depth and lateral coordinates works rather well for sufficiently small strains that cause phase variations of only 1-3 periods over the scan depth (as in Fig. 2 in which 1.5 periods of phase-variation occur along the depth). However, for bigger strains, the quality of such phase compensation noticeably degrades as illustrated in Fig. 4 in which six periods of strain-induced phase-variation are shown. The compensation quality worsening in comparison with the example in Fig. 2 can be explained in the following way. Small strain leads to small displacements of scatterers and most of them do not leave their pixels. Therefore, we can neglect changes in amplitude distribution and the main effect of bulk motion is the small phase variation between B-scans. With increasing strain, more and more scatterers leave their pixels and at the same time the signal from scatterer(s) in a particular pixel gradually leaks to neighboring pixels. This leads to high decorrelation noise in which both phase and amplitude contributions become important.
In view of this, to compensate bulk motions in case of larger strains we developed a new method which uses not only phase compensation but rather combined amplitudephase compensation allowing for a significant reduction of interframe decorrelation. This method is based on the Fourier-shift theorem and is described in the following section in detail.
4 Bulk Motion Compensation Based of Fourier-Shift Procedures
For understanding the principle of the proposed method let us consider a single A-scan first. Amplitude-phase compensation is based on the Fourier shift theorem:
FFT(A(z - Az)) = FFT(A(z))e-kAz. (6)
Therefore, if we need to translationally shift an A-scan, it can be performed as:
Ashiyl (z -Az) = IFFT (FFT (A( z))eikA), (7)
where k is a wavenumber, Az is a required shift in an axial direction, FFT and IFFT - forward and inverse Fourier transforms.
It should be emphasized that it is possible to shift the scan patterns with a subpixel accuracy using Eq. (7). If Az is a multiple of a pixel size Hpx, then the shift is ideal (amplitude and phase do not change their values, just circularly shift) and if Az is not a multiple of Hpx, the complex amplitude somewhat changes as shown in the
right column of Fig. 5. It is important to emphasize that the profile changes not because Fourier shift performs not good enough for subpixel shift but due to the properties of pixelated image formation. Indeed, let us consider a situation when scatterers themselves are translationally shifted by a fraction of a pixel. In this case, point spread function of each scatterer will be sampled at different points than before the shift, which leads to profile changes. Exactly the same situation occurs when we perform subpixel Fourier shift. This effect is illustrated in Fig. 5.
Fig. 4 Simulated example of the phase difference between two subsequent B-scans in case of fairly big strain resulting in the six-period phase variation in the bottom region. In contrast to Fig. 2 where blood vessels are clearly visible, here big strain leads to high decorrelation noise, which masks the signal variability in blood vessel cross-sections. The strain value is 3*10-3. The scale bar is 0.2 mm in both panels.
Az = 2 Hpx Az = 1.5 Hpx
Amplitude
Fig. 5 Amplitude and phase of an A-scan before and after translational shift for different Az. For a shift by an integer number of pixel sizes (left column) the amplitude and phase profiles are simply displaced without changing the form; in contrast, a translational shift by a fractional number of pixels (right column) causes some changes in the phase and amplitude profiles, in other words, the initial and shifted profiles experience decorrelation.
It should be emphasized that it is reasonable to apply the Fourier-shift theorem to the entire scan depth only in case when the scan shift is depth-independent. In case of depth-dependent shift of scatterers, A-scan can be divided into small parts, where the current shift Az is
fairly constant, and therefore, the Fourier shift theorem can be used with a reasonable accuracy. This idea is shown in Fig. 6. It should be emphasized that Eq. (7) leads to a circular shift of the entire array to which it is applied. Therefore, only the central part of the chosen Fourier-shifted fragment of the A-scan is not affected by this effect and can be used for reconstruction of the final shifted image.
Fig. 6 Demonstration of the idea for use of Fourier shift theorem when required shift depends on the depth Az=Az(z). A-scan is divided into parts where Az ~ const and Fourier shift is applicable. The final image is composed of these partial A-scans to which a separate shift was applied.
To demonstrate the method possibilities, two A-scans before and after deformation with constant strain along all depth were generated (see Fig. 7a). In this example the uniform strain is 5 = 3*10-3 and the displacement is zero at the interface of the OCT probe and tissue at z = 0. The required shift Дz for depth z evidently is Д = z ■ s . For such strain, the displacements of scatterers are essentially supra-wavelength within a significant part of the scan and near the bottom have a size of a pixel. Also change in the amplitude profile for the deformed A-scan is also visible. Fig. 7b shows the result of the application of the above-described Fourier-shift procedure with the sliding window of \5 pixels in size shifted pixel by pixel along the deformed A-scan. After this procedure both profiles look much more alike but still slightly differ. The
difference is caused by the fact that when shift is depth-dependent, the used Az is only an approximation and shift cannot be performed ideally. Plus, in the real experiment difference will be somewhat higher because the required shift is not known and has to be estimated using the experimental OCT data. However, even though the difference is still present, the used Fourier-shift procedures significantly improve the quality of OCTA visualization as will be demonstrated in the next Figures.
Now, when the principles of the depth-dependent Fourier-shift procedures are presented, the full sequence of steps of depth-dependent back-shifts to compensate deformation-induced bulk motions can be described as follows.
1. The first step is the initial estimate of the smoothed
A
phase difference A z, x) which coincides with the
A
phase-variation A 9( z, x) defined by Eq. (3). The vectorA
averaged phase difference A 9.nit (z, x) is found using a sliding window. The sizes of the processing window (given by quantities Zp and Xp ) used here should be
smaller than the characteristic scale of deformations but larger than the characteristic scales of the point spread function. Due to such averaging the fairly large-scale phase variations caused by translational and strain-produced displacements can be well estimated and compensated, whereas smaller-scale variations in vessel cross sections will be retained after the compensation.
For illustration, we also calculate the unaveraged phase difference:
АФ»Й (^ x) = arg (B1 (^ (^ x) ) ,
(8)
and compare it with the averaged A qinit (z, x) for the same simulated scans as in Fig. 4. Both unaveraged and averaged phase differences are shown in Fig. 8.
Fig. 7 Compensation of deformation-induced depth-dependent displacements using the Fourier-shift procedures. Panel (a) shows amplitudes of the reference (black) and deformed (red) A-scans and (b) are the reference and deformed A-scans after depth-dependent back-shift. The zoomed insets show the difference/similarity of the compared A-scans in more detail. It is clear that after the back-shift, the two profiles differ much less.
A
Fig. 8 Visualization of phase difference between consecutive B-scans. Panel (a) is pixel-to-pixel unaveraged phase difference A9init (z, x) and (b) is phase
A
difference A qinU(z, x) vector-averaged with a sliding window of 4 x4 pixels in size. Comparison of (a) and (b) demonstrates that the decorrelation noise in the phasevariation distribution is well suppressed, whereas the larger scale strain-induced phase variations are well captured. The scale bar is 0.2 mm in both panels.
2. Based on the calculated smoothed phase difference we calculate the depth-dependent axial shift Az (z, x) caused by the deformation. It can be performed in various ways, e.g., in the most simple case it is done by commonly used unwrap procedure:
(9)
Az(z, x) = unwrap I Д фш (z, x) / 2к0 I.
3. Then we perform the described above Fourier shift using the calculated Az for each A-scan separately using Eq. (7), so that the back-shifted scan B2 shift is obtained
as a result. It should be emphasized that this shift leads only to a shift of the image, but not scatterers themselves. This means that although the image of scatterers is shifted to its "original" position before deformation, its phase corresponds to the deformed position. That means that althugh we moved the complex-valued image to its place before deformation, we still have to perform phase compensation.
4. The next step is calculation of the vector-averaged phase difference once more, but this time between the reference scan B and back-shifted deformed scan
B.
2_sh ift ■
A ф fm ( z, x) = arg
( z+Zp-1 x+Xp-1 Д
X X Bi( z ' x ')
v B*_sMft (z ', x')
(10)
Again, just for comparison we also calculate the unaveraged phase difference:
A9 fin (z, x) = arg (B1 (z, x)K shift (z, x) ). (11)
It is important to emphasize that the above described procedures based on the Fourier shift lead to a significant decrease in the decorrelation noise and, consequently, to smaller amplitude-phase differences between consecutive scans. This effect is shown in Fig. 9.
Fig. 9 Comparison of the phase differences before and after Fourier back-shift for the same simulated scans as in Figs. 4 and 8. Panel (a) is pixel-to-pixel phase difference before the shift, (b) is pixel-to-pixel phase difference after the shift defined by Eq. (11), and (c) is the phase difference between the reference B1 and back-shifted scan B2 shift found with averaging given by
Eq. (10). Comparison of (a) and (b) clearly shows that the proposed Fourier shifting results in significantly smaller decorrelation noise which is further reduced by the moderate averaging as shown (c). The scale bar is 0.2 mm in all panels.
5. The last step is phase compensation for the deformed and back-shifted B-scan:
Bshift( ^ x) = B2_ shift( z, x) x
X exp [iA ф fin( z, x)I.
(12)
Fig. 10 shows the result of the phase compensation given by Eq. (12) in comparison with the phase compensation given by Eqs. (3) and (4) as shown in Fig. (4) in which the back-shift was not used and appreciable decorrelation noise was retained.
Phase compensation Shift compensation
0.2 mm
(b)
Fig. 10 Comparison of the phase difference between B-scans after compensation: phase (a) and shift (b) compensation methods. The proposed method has a much smaller decorrelation noise. The same pair of simulated scans as in Figs. (4), (8), and (9) is used. The scale bar is 0.2 mm in both panels.
л
5 Results of Angiographic Processing of Simulated 3D Data
For comparison of the bulk motion compensation using the phase-only correction and the proposed method based on Fourier-shift, we simulated a 3D stack of B-scans with a channel of moving scatterers imitating a vessel. This vessel has a diameter of 46 ^m and is oriented along the volume diagonal as schematically shown in Fig. \\. The assumed in the simulations rate of obtaining B-scans is 20 Hz (as for the OCT system used in Ref. [20] with the other typical parameters described in Section 2). Erythrocytes are modeled as discrete scatterers flowing in one direction and simultaneously performing Brownian motion (a detailed description of such simulations can be found in Ref. [24]). The regular velocity of moving scatterers was chosen to be \7.4 ^m/s which is quite realistic for small vessels. The size of the image is 1000*400*270 ^m in z, x and y directions, respectively.
Fig. \\ Schematic image of the simulated 3D volume. The diagonal cylinder represents the blood vessel filled with erythrocytes. The rest of the volume also contains scatterers imitating the "solid" tissue but these scatterers are not shown. They have the same concentration and scattering strength as scatterers inside the vessel.
In contact mode the pressure produced by the OCT probe results in appearance of inter-frame strain. It was simulated as inter-B-scan strains with randomly fluctuating amplitude with a Gaussian distribution with zero average and standard deviation STD = 2*10-3. Phase compensation was performed with a 6x6 sliding window, the size of a filtering window in the slow-scanning direction was chosen equal to 7.
The results obtained with the convolutional filtering procedure proposed in Ref. [20] combined with the above-described phase and shift compensations are shown in Fig. 12. They are represented in the traditionally used en face form as maximum amplitude projection of the high-pass filtered 3D data onto the horizontal plane. It is clearly seen that images using the proposed compensation method have much fewer bulk motion-induced artefacts. One of the main advantages of using simulated data for the development of new signal processing methods is that all the parameters of the simulation are known including the position of the vessel, flow velocity, intensity of bulk motions, etc. Therefore, we can reliably quantitatively compare the two compensation methods in highly controllable conditions.
To quantitatively compare the results of different compensation methods we calculate the ratio of the mean signal inside the vessel to the mean signal outside the vessel:
a =
I
/ Npx
vessel ' vessel
ZS / Npx
background background
(13)
Here, Svessel is the amplitude of each of the pixels belonging to the vessel in the angiographic en face image; similarly Sbackground denoted the amplitude of
each pixel outside the vessel; the amount of pixels inside
and outside the vessel are Aground and Nbackground ,
respectively. The summation of the OCTA signals in Eq. (13) is performed over the total number of pixels in the image. For the en face OCTA images shown in Fig. 12, the so-found contrast is a = 5.32 for phase compensation (Fig. 12c) and a = 16.63 for the image obtained with the Fourier shift (Fig. 12c).
Although in Fig. 12 it was supposed that there are no measurement noises other than decorrelation noise, the used simulation method also allows one to readily imitate the noise of the receiving array by adding complex-valued random values to every pixel. For testing of the performance of the two compensation methods in the presence of reception noise we calculated contrast a after adding various noise levels in the simulated scans. The results presented in Table 1 demonstrate that even for rather low signal-to-noise ratios (SNR), the Fourier-shift compensation enables a notably higher contrast. Bearing in mind that typical OCT devices have the reception noise level above 30 dB, significant image quality improvement can be expected for real data.
Table 1 Dependence of the angiographic-image contrast on SNR.
SNR, dB
no noise
40
35
30
25
20
a, phase compensation
5.32
5.31
5.29
5.23
5.06
4.66
a, shift compensation
16.63
16.42 15.93 14.65 12.00
8.44
Fig. 12 Comparison of OCTA images obtained with high-pass filtering using phase (left part) and Fourier-shift (right part) compensations. Panels (a) and (b) are examples of in-depth images; (c) and (d) are en face projections. The proposed Fourier-shift compensation performs noticeably better: there are much fewer artefacts in the background area. The scale bar is 0.05 mm in all panels.
6 Conclusions
Here, we present a new method for bulk motion compensation to improve the quality of OCT-based angiographic visualization in the presence of fairly intense depth-dependent and laterally inhomogeneous strain-induced masking displacements of scatterers in the "solid" tissue, which may occur in the contact-mode angiographic examinations. Although the phase of the pixels is much more sensitive than amplitude to small displacements of scatterers, for elevated strain level, displacements may become so big that scatterers
experience essentially supra-wavelength displacements and may even leave their pixels. This may appreciably affect amplitude profiles of A-scans and lead to high decorrelation noise comprising both phase and amplitude contributions. We present a method intended to significantly reduce such masking strain-induced decorrelation by using the depth-dependent back-shift of scatterers based on the Fourier-shift theorem. Such a compensation method allows one to back-shift an image to its original place before deformation. The proposed procedures lead to lower inter-scan decorrelation, better phase compensation, and therefore, lower signal variability outside blood vessels while preserving high variability inside. A comparison of the earlier used phase-only compensation with the proposed method uses simulated data enabling highly controllable conditions, which is often unachievable in real experiments. The performed comparison shows a much better quality of OCTA image and a significantly higher contrast between the blood vessel and background deformed tissue. It should be emphasized that although the proposed method was tested only on the high-pass filtering method, due to amplitude-phase compensation it can be applicable for any angiographic visualization approach based on amplitude or phase or complex signal variability estimation. Thus, the new masking motion compensation method looks very promising for application to patients, for which utilization of various immobilizing devices is unacceptable and contact-mode OCTA is a realistic acceptable alternative. However, because of residual masking tissue motions obtaining of good-quality angiographic images remains challenging. The proposed improved bulk-motion compensation method will be tested on experimental data in the nearest future.
Disclosures
The authors declare no conflict of interest.
Acknowledgements
The study was supported by the Russian Science Foundation grant No. 22-22-00952.
References
1. V. Y. Zaitsev, A. L. Matveyev, L. A. Matveev, E. V. Gubarkova, A. A. Sovetsky, M. A. Sirotkina, G. V. Gelikonov, E. V. Zagaynova, N. D. Gladkova, and A. Vitkin, "Practical obstacles and their mitigation strategies in compressional optical coherence elastography of biological tissues," Journal of Innovative Optical Health Science 10(06), 1742006 (2017).
2. J.-H. Kang, S.-W. Lee, J.-Y. Yoo, M.-S. Kang, B.-P. Kim, B. S. Yoon, Y. T. Kim, and N.-H. Jo, "Early-Stage Diagnosis of Cervix using the Polarization Sensitive Optical Coherence Tomography(A Preliminary Study)," Korean Journal of Optics and Photonics 18(5), 367-372 (2007).
3. A. Weatherbee, M. Sugita, K. Bizheva, I. Popov, and A. Vitkin, "Probability density function formalism for optical coherence tomography signal analysis: a controlled phantom study," Optics Letters 41(12), 2727 (2016).
4. K. V. Larin, D. D. Sampson, "Optical coherence elastography-OCT at work in tissue biomechanics," Biomedical Optics Express 8(2), 1172-1202 (2017).
5. 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).
6. A. A. Plekhanov, M. A. Sirotkina, A. A. Sovetsky, E. V. Gubarkova, S. S. Kuznetsov, A. L. Matveyev, L. A. Matveev, E. V. Zagaynova, N. D. Gladkova, and V. Y. Zaitsev, "Histological validation of in vivo assessment of cancer tissue inhomogeneity and automated morphological segmentation enabled by Optical Coherence Elastography," Scientific Reports 10(1), 11781 (2020).
7. E. V. Gubarkova, A. A. Sovetsky, V. Yu. Zaitsev, A. L. Matveyev, D. A. Vorontsov, M. A. Sirotkina, L. A. Matveev, A. A. Plekhanov, N. P. Pavlova, S. S. Kuznetsov, A. Yu. Vorontsov, E. V. Zagaynova, and N. D. Gladkova, "OCT-elastography-based optical biopsy for breast cancer delineation and express assessment of morphological/molecular subtypes," Biomedical Optics Express 10(5), 2244 (2019).
8. A. Mariampillai, B. A. Standish, E. H. Moriyama, M. Khurana, N. R. Munce, M. K. K. Leung, J. Jiang, A. Cable, B. C. Wilson, I. A. Vitkin, and V. X. D. Yang, "Speckle variance detection of microvasculature using swept-source optical coherence tomography," Optics Letters 33(13), 1530 (2008).
9. E. Jonathan, J. Enfield, and M. J. Leahy, "Correlation mapping method for generating microcirculation morphology from optical coherence tomography (OCT) intensity images," Journal of Biophotonics 4(9), 583-587 (2011).
10. C.-L. Chen, R. K. Wang, "Optical coherence tomography based angiography [Invited]," Biomedical Optics Express 8(2), 1056 (2017).
11. T. E. de Carlo, A. Romano, N. K. Waheed, and J. S. Duker, "A review of optical coherence tomography angiography (OCTA)," International Journal of Retina and Vitreous 1(1), 5 (2015).
12. U. A. T. Hofmann, J. Rebling, H. Estrada, P. Subochev, and D. Razansky, "Rapid functional optoacoustic microangiography in a burst mode," Optics Letters 45(9), 2522 (2020).
13. V. X. D. Yang, M. L. Gordon, B. Qi, J. Pekar, S. Lo, E. Seng-Yue, A. Mok, B. C. Wilson, and I. A. Vitkin, "High speed, wide velocity dynamic range Doppler optical coherence tomography (Part I): System design, signal processing, and performance," Optics Express 11(7), 794 (2003).
14. R. A. Leitgeb, L. Schmetterer, W. Drexler, A. F. Fercher, R. J. Zawadzki, and T. Bajraszewski, "Real-time assessment of retinal blood flow with ultrafast acquisition by color Doppler Fourier domain optical coherence tomography," Optics Express 11(23), 3116 (2003).
15. R. K. Wang, L. An, "Doppler optical micro-angiography for volumetric imaging of vascular perfusion in vivo," Optics Express 17(11), 8926 (2009).
16. D. M. Schwartz, J. Fingler, D. Y. Kim, R. J. Zawadzki, L. S. Morse, S. S. Park, S. E. Fraser, and J. S. Werner, "Phase-Variance Optical Coherence Tomography: A Technique for Noninvasive Angiography," Ophthalmology 121(1), 180-187 (2014).
17. S. Makita, Y. Hong, M. Yamanari, T. Yatagai, and Y. Yasuno, "Optical coherence angiography," Optics Express 14(17), 7821 (2006).
18. Y. Jia, O. Tan, J. Tokayer, B. Potsaid, Y. Wang, J. J. Liu, M. F. Kraus, H. Subhash, J. G. Fujimoto, J. Hornegger, and D. Huang, "Split-spectrum amplitude-decorrelation angiography with optical coherence tomography," Optics Express 20(4), 4710 (2012).
19. J. Tang, S. E. Erdener, S. Sunil, and D. A. Boas, "Normalized field autocorrelation function-based optical coherence tomography three-dimensional angiography," Journal of Biomedical Optics 24(03), 1 (2019).
20. A. Moiseev, S. Ksenofontov, M. Sirotkina, E. Kiseleva, M. Gorozhantseva, N. Shakhova, L. Matveev, V. Zaitsev, A. Matveyev, E. Zagaynova, V. Gelikonov, N. Gladkova, A. Vitkin, and G. Gelikonov, "Optical coherence tomography-based angiography device with real-time angiography B-scans visualization and hand-held probe for everyday clinical use," Journal of Biophotonics 11(10), e201700292 (2018).
21. V. Demidov, L. A. Matveev, O. Demidova, A. L. Matveyev, V. Y. Zaitsev, C. Flueraru, and I. A. Vitkin, "Analysis of low-scattering regions in optical coherence tomography: applications to neurography and lymphangiography," Biomedical Optics Express 10(8), 4207 (2019).
22. V. Yu. Zaitsev, I. A. Vitkin, L. A. Matveev, V. M. Gelikonov, A. L. Matveyev, and G. V. Gelikonov, "Recent Trends in Multimodal Optical Coherence Tomography. II. The Correlation-Stability Approach in OCT Elastography and Methods for Visualization of Microcirculation," Radiophysics and Quantum Electronics 57(3), 210-225 (2014).
23. L. A. Matveev, V. Yu. Zaitsev, G. V. Gelikonov, A. L. Matveyev, A. A. Moiseev, S. Yu. Ksenofontov, V. M. Gelikonov, M. A. Sirotkina, N. D. Gladkova, V. Demidov, and A. Vitkin, "Hybrid M-mode-like OCT imaging of three-dimensional microvasculature in vivo using reference-free processing of complex valued B-scans," Optics Letters 40(7), 1472 (2015).
24. 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).
25. G. Yao, L. V. Wang, "Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media," Physics in Medicine & Biology 44(9), 2307-2320 (1999).
26. M. Kirillin, I. Meglinski, V. Kuzmin, E. Sergeeva, and R. Myllyla, "Simulation of optical coherence tomography images by Monte Carlo modeling based on polarization vector approach," Optics Express \8(2\), 2\7\4 (20\0).
27. A. Doronin, I. Meglinski, "Online object oriented Monte Carlo computational tool for the needs of biomedical optics," Biomedical Optics Express 2(9), 246\ (20\\).
28. 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 \\(\0), \0560\ (20\4).
29. A. L. Matveyev, L. A. Matveev, A. A. Moiseev, A. A. Sovetsky, G. V. Gelikonov, and V. Y. Zaitsev, "Semi-analytical full-wave model for simulations of scans in optical coherence tomography with accounting for beam focusing and the motion of scatterers," Laser Physics Letters \6(8), 085601 (20\9).
30. M. Almasian, T. G. van Leeuwen, and D. J. Faber, "OCT Amplitude and Speckle Statistics of Discrete Random Media," Scientific Reports 7(\), \4873 (20\7).
