1.

Introduction

Optical coherence tomography (OCT), as a noninvasive diagnostic tool, has been widely used in dermatological research and for clinical skin diagnoses for decades now.^1–5 In addition to examining skin morphology, important advances have been made toward the imaging of functional structures, such as the vascular networks, within a tissue depth of $\sim 1\mathrm{mm}$. Such developments have been achieved using a process known as optical coherence tomography angiography (OCTA). OCTA utilizes an endogenous contrast mechanism that originates from the fluctuations of blood flow signals to differentiate a vascular network from its encapsulating static tissue bed.^6–9 To extract the signal fluctuations, B-scans are typically repeated at one cross-sectional position to create a cluster of comparable B-frames.^6,10 Signals coming from static tissue within such B-frames show little interframe variance, while those from blood vessels experience severe fluctuations due to the continuous flow of blood cells. To aid in visualization, three-dimensional (3-D) vascular images are commonly constructed from a series of blood flow cross sections by raster scanning across the slow transverse direction. En face view vasculatures can be visualized through the mean intensity or maximum intensity projection of blood flow signals at a given depth.

Practically speaking, the process of 3-D OCTA imaging of human skin in situ is very susceptible to motion-derived artifacts stemming from various sources, such as respiratory or cardiac activities, involuntary contracting by the subject, and unintentional movements of the handheld probes during dermatological scanning. From previous reports,^11–14 motion-derived artifacts resulted in the exclusion of sizable volumes of data from analyses. Even relatively minor movements, such as those attributed to the heartbeat, microtremors in the skin, or motions due to breathing, can increase background noise, predominantly resulting in the appearance of line artifacts in the en face projection. Fortunately, coregistering the B-frame clusters, i.e., repeated B-scans at the same spatial location, before applying angiography algorithms can help with the reduction of such artifacts. Within individual B-frames, the current OCT speeds, typically with an A-line rate of $\sim 100\mathrm{kHz}$ and a B-frame rate of $\sim 180\mathrm{fps}$, usually limit the level of artifacts seen from severe movements, e.g., inadvertent sliding of the probe. However, the same cannot be said for consecutive frames because the time between adjacent frames is relatively long.

Reported efforts on compensating for motion-derived artifacts in OCTA can be summed up within two categories: the hardware approach and the software approach. Hardware techniques adopting high scan rates, e.g., megahertz OCT,¹⁵ are proposed to reduce the acquisition time during in vivo skin imaging, and consequently the probability of motion involvement is lessened. In commercialized ophthalmology OCTA systems,¹⁶ for example, scanning laser ophthalmoscopy^17,18 is a prevalent strategy used for detecting rapid/large motions and guiding appropriate rescans with specified overlapping to counter them. Adopting such a technique for imaging of the skin, however, would prove fruitless, mainly because unlike retinal imaging, human skin has much higher scattering properties compared to the relatively transparent ocular tissue. Furthermore, the increased complexity and cost associated with acquiring and developing additional devices may reduce its competitiveness in the application of the technology in dermatology.

An alternative solution is software-aided postprocessing through image registration,¹⁹ in which two or more datasets are usually aligned and merged into one so that artifacts can be replaced by their spatial counterparts where there is no apparent artifact. Using this concept, Kraus et al.²⁰ corrected for motion-induced distortion by retrospectively coregistering orthogonally acquired volumetric scans. Although initially encouraging, this approach required considerable computational power to register just two datasets. Another approach attempted by Antony et al.²¹ successfully extracted tissue features from OCT structure images based on surface segmentation, and accordingly, they designed geometric transformations to correct the transverse distortions in 3-D. A fiducial marker specifically designed for human skin imaging has also been utilized to achieve 3-D motion correction on an individual B-frame basis.²² Although this method did establish good image correspondence, the feature-based OCT structure registration finds it difficult to satisfy the accuracy requirement for microscale angiography. Moreover, if directly applying the feature-based method to register blood flow signals, one usually has difficulties in assessing the essential features of OCT angiograms, such as extreme curvatures and/or contours. The fiducial marker adhesively bonded to the skin may also cause an allergic response or inflammation that not only leads to future reluctance from the patients to rescan but also interferes with dermatopathological evaluations. An interesting study carried out by Hendargo et al.²³ first presented an automatic nonrigid registration for use in OCTA by directly manipulating the orthogonal en face retinal vasculatures. Additionally, Zang et al.²⁴ further demonstrated the practicability of nonrigid registration by minimizing stitching artifacts when montaging en face vascular images from equidirectional scanned volumes.

Here, we present a fully automated motion compensation method that employs both rigid transformation and nonrigid free-form deformation to correct motion-derived artifacts that frequently occur in in vivo human skin OCTA imaging. Although orthogonal transverse scans are widely used in ophthalmology research, the restricted uniform sampling therein becomes inappropriate when imaging irregular tissue characteristics, e.g., an incision in the skin. Therefore, a nonorthogonal dual-volume scanning procedure was utilized within our protocol. A rigid affine transformation was first applied to the reticular vasculature of the skin, which primarily consists of deep, large vessels. The generated affine parameters were then applied to the whole vessel network with additional translation adjustment through rigid subpixel compensation. In addition, a B-spline free-form deformation was used to calibrate localized misalignment due to nonrigid distortions and accurately register the dual vasculatures into one artifact-free composition. The whole registration process was feedback controlled by vascular image quality evaluations. The overall performance of the motion compensation algorithm, in terms of correlation coefficient (CC) and peak signal-to-noise ratios (PSNR) of the registered vasculatures, was evaluated and compared to that prior to processing. Our intention is to provide a technique that will offer greater opportunities for a number of clinical outcomes, such as the improved early diagnosis of skin abnormalities or the automatic quantification of dermatosis severity, for example.

2.

Methods

2.1.

Imaging System and Data Acquisition

To scan human skin, we employed an in-house-built swept source OCT (SS-OCT) system. The in vivo imaging study that uses home-built systems to image human subjects was reviewed and approved by the Institutional Review Board of University of Washington, and informed consent was obtained in writing from all subjects before imaging. This study followed the tenets of the Declaration of Helsinki and was conducted in compliance with the Health Insurance Portability and Accountability Act.

Briefly in the SS-OCT, a 100-kHz vertical-cavity surface-emitting swept laser source (SL1310V1-10048, Thorlabs Inc., Newton, New Jersey) was used as the light source providing a central wavelength of 1300 nm and a spectral bandwidth of 68 nm. Light from the laser was first split into sample and reference arms through a $90/10$ coupler. In the sample arm, a $5\times$ objective lens (LSM03, Thorlabs Inc.) focused a beam spot into the skin with an incident power of 5 mW. The beam spot was scanned by a paired $X\text{-}Y$ galvo scanner, forming raster sampling patterns that consisted of a series of fast ($x$-axis) and slow ($y$-axis) scans. The sample arm was configured as a handheld probe affixed with a sample spacer (central hole $Ø=25\mathrm{mm}$) to keep a consistent distance between the objective lens and the skin. The interference signal was detected by a 1.6-GHz photodetector (PDB480C-AC, Thorlabs Inc.) and sampled by a 12-bit digitizer (ATS9350, Alazartech) at $500\mathrm{MS}/\mathrm{s}$. The resulting axial resolution was measured to be $15\mu \mathrm{m}$ in air, and the transverse resolution was measured to be $15\mu \mathrm{m}$ at the focus. The system sensitivity was measured at 105 dB.

Dual-volume datasets were sequentially acquired through a specifically designed C-scan protocol, as shown in Fig. 1(a), in which each volume consisted of a total of 2500 B-scans with five repetitions at one spatial location, leading to 500 B-scan positions. The C-scans were acquired over two volumes, $7\times 7\times 12{\mathrm{mm}}^3$ in each, with a constant sampling voxel size of $14\times 14\times 10\mu {\mathrm{m}}^3$. The total acquisition time for two volumes was $\sim 30\mathrm{s}$.

Fig. 1

C-scan protocol and resultant cross sections: (a) dual-volume OMAG scanning protocol, (b) a typical cross-sectional structure image, and (c) corresponding blood flow image.

Before registering the volume datasets, the acquired frames were preprocessed through the following steps. First, we coregistered the repeated B-scan frames through a subpixel image registration algorithm²⁵ to minimize any small or slow motions. Thereafter, an optical microangiography (OMAG) algorithm⁷ was applied to extract dynamic blood flow information from the static tissue background. A typical cross-sectional structure image and its corresponding blood flow are presented in Figs. 1(b) and 1(c), respectively. Second, an automatic segmentation approach was adopted from Yin et al.²⁶ to differentiate the papillary dermis layer, filled with an abundant superficial plexus, and the reticular dermis layer, which is dominated by larger arterioles or venules, as shown in Figs. 2(a), 2(d), and 2(g), respectively. Furthermore, the en face vasculatures were separated into two images according to the scanning protocol, as shown in Figs. 2(b) and 2(c), 2(e) and 2(f), or 2(h) and 2(i). Severe motion-derived artifacts are visualized in Fig. 2 as vertical lines marked by red arrows.

Fig. 2

En face skin vasculatures from different depths. (a) Motion affected en face arm skin vasculature obtained by two repeated C-scans with (b) and (c) representing single-volume datasets. (d) Dual-scanned en face skin vasculature of papillary layer with (e) and (f) representing single-volume datasets. (g) Dual-scanned en face skin vasculature of reticular layer with (h) and (i) representing single-volume datasets. Red arrows indicate severe motion-derived artifacts.

2.2.

Image Registration Method

During image registration, two different approaches were utilized sequentially based upon the vascular network within the dual-volume datasets: a general registration that was mainly achieved by rigid affine registration plus rigid subpixel motion compensation and a precise registration that was mainly realized through nonrigid B-spline free-form deformation, as shown in Fig. 3 (dashed boxes). This section will outline in detail the step-by-step motion compensation procedure applied for human skin OCTA. The preregistration process introduced in Sec. 2.1 is summarized here as the elliptical blocks of Fig. 3.

Fig. 3

Flow chart outlining automatic motion compensation in skin OCTA. $P$ represents the papillary vasculature, $R$ represents the reticular vasculature, and $\alpha$ represents the rescaling factor. Function of each process step is classified in the legend.

2.2.1.

Peak detection-based frame division

After separation of the dual scans, the en face vasculatures of both depths, the papillary layer and the reticular layer, were successfully acquired, but with severe motion-derived artifacts that appear as white lines. The locations of those artifacts were identified in the all-depth vasculatures by averaging pixel intensities along the $x$-axis, as shown in the one-dimensional (1-D) projection curves in Figs. 4(a) and 4(b). Peak detection according to the peak prominence was used to identify which of the projected B-frames had suffered from motion-derived artifacts. The peak prominence was empirically selected as 0.3 times the maximum peak elevation above the mean intensity (mean refers to the average intensity of the whole image). The calculated peak widths were also utilized to indicate how many frames were affected by each motion, which were then accordingly removed. When two adjacent peaks were detected with overlapping widths, as shown in the dashed box of Fig. 4(a), they were automatically regarded as one peak to avoid confusion in the following frame division. Each vasculature was then divided at the location of the peaks, generating a sequence of frame groups that were free of artifact lines. Last, the peak detection results were extended to both the papillary layer and the reticular layer for registrations in the next step. This approach has several advantages over the previously published threshold-based motion detection,^23,24 especially for human skin OCTA. As illustrated in Fig. 4(e), when local vascularity increases to a level beyond that of normal tissue, e.g., angiogenesis during wound healing, adopting a simple threshold technique could deliver inaccurate results by either mistakenly regarding a normal region as an artifact or by erroneously taking an artifact as true data. Moreover, the peak width also provided a better indicator for the occupied region of each specific motion-derived artifact compared with that using arbitrary thresholds.

Fig. 4

Peak detection-based frame division. (a) 1-D projection curve showing artifact peaks in vasculature $A$ with overlapping peak widths marked in a dashed box. (b) 1-D projection curve showing artifact peaks in vasculature $B$. Detected peaks are marked with black dots. (c) Vasculature $A$ with artifact line removal and frame division. (d) Vasculature $B$ with artifact line removal and frame division. (e) An example showing the advantages of peak detection over threshold in detecting motion-derived artifacts. The threshold approach fails by regarding normal regions as artifacts (upper subgraph, threshold 1) or regarding artifact lines as normal (lower subgraph, threshold 2) due to regional signal variance.

2.2.2.

Sequential general registration of frame groups

The divided groups of B-frames within both reticular layer volume scans were labeled in ascending order according to their original physical positions. Both reticular volumes were then sequentially coregistered into one through two-dimensional (2-D) affine transformation,²⁷ as shown in Fig. 5(a). To represent the affine transformation with matrix operations, we adopted homogeneous coordinates in which 2-D coordinates along the $x$- and $y$-axes $[x,y]$ were represented by 3-D coordinates $[x,y,1]$. Therefore, typical motion-induced distortions, such as shift, expand/contract, twist, or rotation, within the frame groups can be corrected through the combination of translation, scale, shear, or rotation transformations, as described in Figs. 5(b)–5(e), respectively.

Fig. 5

Affine registration of the reticular vasculature. (a) Schematic registration procedure of reticular vasculature. A_R represents the first reticular vasculature, B_R represents the second reticular vasculature, and Reg_R represents the affine registered reticular vasculature. (b) to (e) Affine transformation diagrams and matrix operations accounting for translation, scale, shear, and rotation, respectively. $S$ represents the reference vasculature as registration model and $M$ represents the moving vasculature to be registered.

The affine registration was evaluated by a gradient descent method, the limited memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS), which maximized the normalized mutual information between en face vasculatures of registered frame groups^27,28

Eq. (1)

$${\text{Cost}}_{MI}=\frac{H(S)+H[T_r(M)]}{H[S,T_r(M)]},$$

where $H(S)$ represents the marginal entropy of static vasculature $S$ as the registration model, $H[T_r(M)]$ represents the marginal entropy of moving vasculature $M$ to be registered, and $H[S,T_r(M)]$ represents their joint entropy, as defined by Maes et al.²⁹ The cost function reaches maximum when the two vasculatures are well aligned.

During the affine registration of the reticular layer, affine parameters for each frame group were recorded and thereafter applied to the en face vasculatures of both depths, as expressed in Fig. 6(a). After each affine transformation, a rigid subpixel registration algorithm was added to compute the upsampled ($10\times$) cross correlation between the images to be registered and the reference static images by means of fast Fourier transform and correlation peak positioning; accordingly, subpixel shifts down to $1/10$-pixel scale were compensated for. The computational burden for this subpixel compensation was acceptable when considering the sample size of $500\times 500$ for each vasculature. Notice that the all-depth vasculatures used here in general registration were rescaled as combinations of the papillary projection plus a visualization enhanced reticular projection, as shown in Fig. 6(b). The enhancement factor, $\alpha$, was empirically established as 2.0 to boost the weight of the large vessels in the registration in comparison to that of the small vessels. The generally registered vasculature composite, Gel_Reg, combining all the transformed frame groups, is shown in Fig. 6(c). Prior to precise registration, it was necessary to separately stitch the transformed frames (from vasculature $A$ or $B$, respectively) into two vasculatures with the motion-derived artifact lines replaced by their counterparts in Gel_Reg. With that, two artifact-reduced images were obtained, denoted as Gel_Reg_A and Gel_Reg_B, as displayed in Fig. 6(c).

Fig. 6

Sequential general registration of the all-depth vasculature. (a) Transform equation to describe the application of affine parameters. (b) Intensity rescale of vasculature $A$ and $B$. AD represents the all-depth vasculature, $P$ represents the papillary vasculature, $R$ represents the reticular vasculature, and $\alpha$ represents the scale factor used to boost visualization of the reticular vasculature. (c) Schematic diagram of the general registration process. Gel_Reg represents the generally registered vasculature, Gel_Reg_A and Gel_Reg_B represent the generally registered vasculatures $A$ and $B$ with artifact lines replaced by counterparts in Gel_Reg.

2.2.3.

B-spline free-form transformation for precise registration

Following sequential general registration, the two en face vasculatures were ready for coregistration. Skin is naturally highly susceptible to significant subject-dependent variations and microvasculature deformations. With that, such deformations cannot be described as parameterized transformations, and so, to compensate for this, we chose to apply a registration algorithm that uses nonrigid cubic B-spline free-form transformation.^30–32

As shown in Fig. 7(a), a typical 1-D cubic B-spline is made up of four individual parts and in each, the curve can be expressed as a third-order polynomial function, as follows:

Eq. (2)

$$\begin{gathered} B_2(t_2)=t_2^3/6 \\ {B}_1(t_1)=(-3t_1^3+3t_1^2+3t_1+1)/6 \\ {B}_0(t_0)=(3t_0^3-6t_0^2+4)/6 \\ {B}_{-1}(t_{-1})={(1-t_{-1})}^3/6. \end{gathered}$$

The parameters, $t_2$, $t_1$, $t_0$, and $t_{-1}$, vary from 0 (left) to 1 (right), and therefore, can be written as $t_2=\frac{x+2u}{u}$, $t_1=\frac{x+u}{u}$, $t_0=\frac{x}{u}$, $t_{-1}=\frac{x-u}{u}$, with $u$ representing the part length, i.e., the distance between adjacent control points located at $-64$, $-32$, 0, 32, and 64. This B-spline curve has two main properties: first, it is continuous at the connecting points between parts and at these points and has continuous first- and second-order derivatives. This characteristic provides superior smoothness for the B-spline curves and for any combinations or interpolations that are based on similar B-splines. For instance, in Fig. 7(b), the value at a given position (green arrow) can be interpolated by four neighboring B-splines curves, which are scaled by the corresponding control points (dashed arrows), as expressed by

Eq. (3)

$$T(x,y)=\sum _{i=-1}^2{B}_i(t_i){\mathrm{\Psi }}_i,$$

where ${\mathrm{\Psi }}_i$ represents the values of control points. Second, a single B-spline curve is nonzero only in a small region; for example, the B-spline in Fig. 7(a) has nonzero values only from $-64$ to 64. This means that any changes based on B-spline are highly localized, so only limited surrounding points are influenced by a single control point, and the computational efficacy of B-spline-based interpolation is acceptable, even for a large number of control points. By comparing the dash curves in Figs. 7(b), $C_1$, and 7(c), $C_2$, the decreased value at control point 0 only leads to a decreased distribution around 0 without extending to regions under $-64$ or beyond 64.

Fig. 7

Schematic explanation of cubic B-spline transformation. (a) A typical 1-D cubic B-spline with four pieces (piece length $\mu$ equals to 32) described by third-order polynomial functions $B_2(t_2)$, $B_1(t_1)$, $B_0(t_0)$, $B_{-1}(t_{-1})$. (b) Schematic diagram of B-spline interpolation with a green arrow denoting where a value to be interpolated was scaled by the dashed arrows denoting four sounding control points. The transformed 1-D B-spline is visualized as a dash curve ($C_1$). (c) Schematic diagram to show the localization of B-spline transformation with the sudden change of a control point at 0. The transformed 1-D B-spline is visualized as a dash curve ($C_2$). (d) Schematic diagram of 2-D B-spline free-form deformation.

Further to that, 2-D B-splines can be simply derived from the products of two orthogonal 1-D B-splines and still possess the characteristics of smoothness and local controllability. Instead of using these weighted summations to interpolate pixel intensity, however, they could be used to characterize pixel displacement from the original position to the warped point; thereby, deforming the distribution of surrounding pixels through B-spline transformation. The transformation equation for each pixel can be described as a linear superposition of 2-D B-splines scaled by control points ${\mathrm{\Psi }}_{h,c}$, which have two subscripts, $h$ and $c$, indicating their positions within a grid of $H_x\times C_y$ pixels

Eq. (4)

$$T(x,y)=\sum _{j=-1}^2\sum _{i=-1}^2{B}_j(t_j)B_i(t_i){\mathrm{\Psi }}_{h+i,c+j},$$

in which $i$, $j$, $t_i$, and $t_j$ are represented by transverse coordinates $x$, $y$ as $h=\lfloor x/H_x\rfloor -1$, $c=\lfloor y/C_y\rfloor -1$, $t_i=\lfloor x/H_x\rfloor -\lfloor x/H_x\rfloor$, $t_j=\lfloor y/C_y\rfloor -\lfloor y/C_y\rfloor$. $\lfloor \rfloor$ is a flooring operation toward the nearest integer. Consequently, the basic concept of 2-D B-spline registration comes down to designing a grid of control points with third-order basic splines, and using the combinations of those splines to adjust misaligned images. As schematically described in Fig. 7(d), the original image depicted by the black mesh grid can be deformed to the blue mesh grid, point by point, according to the transformation operator $T(x,y)$.

After general registration, the vasculatures, Gel_Reg_A and Gel_Reg_B, were zero-padded all around with an extra 50 pixels (each vasculature reaching a sample size of $600\times 600$). Without loss of generality, Gel_Reg_B, being the moving image, was subjected to 2-D B-spline transformation to coregister with Gel_Reg_A, the static reference image. In the interest of calculation efficacy, a multiresolution gradient descent approach was utilized through the L-BFGS optimizer to refine the registration, starting from a grid of $8\times 8$ control points to a grid of $22\times 22$ points. To constrain the vasculature smoothness, a penalty term was applied to regulate the transformation function, expressed as

Eq. (5)

$$\text{Regulation}=\frac{1}{XY}\underset{y=0}{\overset{Y}{\int }}\underset{x=0}{\overset{X}{\int }}\{{\left[\frac{{\partial }^{2}T(B)}{\partial x^2}\right]}^2+{\left[\frac{{\partial }^{2}T(B)}{\partial y^2}\right]}^2+{\left[\frac{{\partial }^{2}T(B)}{\partial xy}\right]}^2\}\mathrm{d}x\mathrm{d}y,$$

where $T(B)$ represents the transformed moving image, and $X$ and $Y$ are the evaluated pixel numbers along the $x$- and $y$-axes, respectively. Additionally, a similarity criterion was designed to quantify the degree of alignment between the static image and the moving image, as expressed in the following sum of squared differences³³

Eq. (6)

$$\text{Cost}=\frac{1}{XY}\underset{y=0}{\overset{Y}{\int }}\underset{x=0}{\overset{X}{\int }}{[A-T(B)]}^2\mathrm{d}x\mathrm{d}y,$$

where $A$ represents the static vasculature.

The B-spline coefficients were then searched out by repeatedly minimizing a combined evaluation function (EF) as

Eq. (7)

$$EF=-\text{Cost}+\beta ·\text{Regulation},$$

where $\beta$ works as a weighting factor to balance the registration accuracy and vasculature smoothness, which was rationally chosen as $\beta =0.01$, based on the size of the current control mesh grids. If no desirable registration coefficient was identified, the searching would automatically stop by exceeding the iteration number and reaching the maximum permissible grid size. In that case, the best approximations would be returned as the proper registration parameters. Here, a composite image was created by summing both the static image and the transformed moving image with equal weights. The composite image was further trimmed by removing previously zero-padded regions.

2.2.4.

Precise registration feedback controlled by vasculature evaluation

The effectiveness of the general and precise registration protocols was evaluated by a specifically designed vasculature evaluation process, as shown in Fig. 8. Taking the generally registered images, Gel_Reg_A and Gel_Reg_B, as an example, two images were normalized to the same mean intensity for objective comparisons to be made. A vessel subtraction map (VSM) was then obtained by subtracting Gel_Reg_B from Gel_Reg_A, as in Fig. 8(a). The VSM was subsequently segmented into $10\times 10$ blocks with the standard deviations between the vessels being calculated from pixel value variations within each block. Finally, a vasculature evaluation map (VEM) was obtained through bicubic interpolation of the blocks to form a regional evaluation criterion with the same sampling as the original vasculatures ($500\times 500$).

Fig. 8

Vasculature evaluation process to assess the quality of registrations. (a) VSM obtained through the normalized subtraction between two coregistered vasculatures. Gel_Reg_A represents the generally registered vasculature $A$ and Gel_Reg_B represents the generally registered vasculature $B$. (b) VEM obtained by calculating regional standard deviation of the VSM and smoothing through a bicubic interpolation. $(n,m)$ represents a typical pixel in the VSM, $N$ and $M$ represent the evaluated pixel numbers along $x$- and $y$-axes.

The VEM was calculated for the generally registered vasculatures, $A$ and $B$, as shown in Fig. 9(a), with the associated vasculature composite ($A+B$), shown in Fig. 9(b). Subsequently, a second VEM was obtained by comparing the B-spline deformed vasculature $B$ from the first-round precise registration with vasculature $A$, as shown in Fig. 9(c), as well as the corresponding composite image, as shown in Fig. 9(d). A significant improvement in image quality was validated by the evaluation maps (see the reduced color bar values) and coregistered vasculatures [see the comparative vessel patterns indicated by the red arrows in Figs. 9(b) and 9(d)]. The precise registration, however, was not conclusively validated due to the relatively coarse B-spline interpolation grid (maximum $22\times 22$ points), which resulted from a trade-off with calculation efficiency. In its place, we used the VEM and identified the regions with mean values exceeding 50% of the maximum VEM value. Those regions within both vasculatures were then segmented out and zero-padded to their closest rectangle blocks. Then, a bicubic interpolation was applied to upsample those blocks twice along both $x$- and $y$-axes. Finally, the interpolated images were rerun through B-spline registration, and afterward downsampled back to the original transverse spacing. This feedback controlling process was repeated several times until reaching the limits where no such regions were detected, as illustrated in the precise registration element of Fig. 3. Using this method, the VEM values across the whole map were $<0.04$, as shown in Fig. 9(e). The corresponding final composite vasculature is shown in Fig. 9(f).

Fig. 9

Vasculature evaluation maps and associated composite vasculatures postregistration. (a), (c), and (e) Vasculature evaluation maps for the general registration, the first-round precise registration, and the second-round precise registration, respectively. (b), (d), and (f) The composite vasculatures obtained after the general registration, the first-round precise registration, and the second-round precise registration, respectively.

3.

Results

Visually, motion compensation and artifact reduction were successfully achieved through step-by-step rigid affine registration, rigid subpixel registration, and nonrigid B-spline free-form registration, feedback controlled by vasculature evaluation. Additionally, a considerable signal-to-noise improvement can be seen when comparing the registered vasculature composite, shown in Fig. 10(b), with the two original vasculatures, shown in Fig. 10(a). Two additional cases are shown in Figs. 10(d) versus 10(c) and Figs. 10(f) versus 10(e). As much data as possible from each vascular scan were carried forward to the artifact-free composite, as demonstrated in Figs. 10(a)–10(f) by arrows. To quantitatively evaluate the results of the registration process and to assess how much the motion compensation algorithm improved the quality of the vasculature, the CC and PSNR of both vasculatures, pre- and postregistrations, were calculated using Eqs. (8) and (9)³⁴

Fig. 10

Representative skin vasculatures before and after motion compensation. (a) Two en face vasculatures of case 1 preregistration. (b) Final artifact-free vasculature of case 1 postregistration. (c) Two en face vasculatures of case 2 preregistration. (d) Final artifact-free vasculature of case 2 postregistration. (e) Two en face vasculatures of case 3 preregistration. (f) Final artifact-free vasculature of case 3 postregistration. Red arrows highlight artifact lines in vasculatures preregistration. Zoomed images within boxes highlight positions with signal-to-noise improvement.

Table 1

CC and PSNR values of vasculatures pre- and postregistration.

4.

Discussion

The proposed registration method was successfully utilized to eliminate motion-derived artifacts in human skin OCTA. While it is common practice to use an edge enhancement filter, such as Gabor filter or Hessian filter, to increase visual contrast (common approaches for visual enhancement³⁶) prior to registration, no such filter was applied to the angiograms in this study. The decision was made to forego such a routine because it could have resulted in the misrepresentation of vessel textures. For instance, the Gabor filter is known to cause extensions of the vessel when used to smooth the vessel edge, which inevitably raises the quantification values of vessel area density and vessel flux. Such an occurrence can happen when capillaries or small-scale vessels are heavily involved in an evaluation, such as when imaging patients with diabetic dermopathy, necrobiosis lipoidica diabeticorum, or pruritus. Possible misreading to any degree would have been unacceptable here as accuracy is the key to this study.

A more critical step in the motion compensation process, however, is accurately identifying the motion-derived artifacts. As discussed in Sec. 2.2.1, peak detection for frame group division is considered a more accurate technique when compared to arbitrary thresholds. This is justified if one considers peak analysis as peak prominence and peak width performing like dynamic thresholds on specific artifacts with consideration given to regional features. Additionally, such a technique allows for the detection of incomplete or mild artifacts. Moreover, to ensure the robustness of this technique, when no artifact was found to exceed our default prominence level, the highest peak was regarded as the point where frame groups would be divided. Although the default prominence of 0.3 times maximum elevation worked well for all cases here, further optimization could be achieved with the use of more regional characteristics, if desired. It should be noted this technique is not without a potential flaw. A conceivable problem exists in mistaking several frames as motion-derived artifacts when large vessels are orientated parallel to those frames, leading to high OMAG signals. Fortunately, however, this problem very rarely occurred here, for two reasons: first, the vessels present in the skin comprise high network distributions that typically would not be captured in single frames, and second, large unidirectional vessels in the skin typically reside in the hypodermis layer, which is almost inaccessible to the OCT beam.

During affine registration, corresponding pixels were measured for statistical dependence via normalized mutual information, which is assumed to be maximal if the images are geometrically aligned. The mutual information equation was derived from the information theory, which calculates the amount of information that one image, $I_1$, contains about another, $I_2$, as³⁷

Compared with the sum of squared differences used in the precise registration, the mutual information is proposed to be less sensitive to image intensity, and therefore, provided a certain robustness for the affine registration. This is especially important here because we did not apply histogram equalization to the original vasculatures, as its nonlinear manipulation may bias the quantitative analyses that follow.³⁸ Additionally, to avoid any dependency of the registration on the vasculature overlap, an improved normalized mutual information²⁸ was used as the standard for image alignment, as in Eq. (1). While, in practice, the affine registration still requires the two vasculatures to preserve some degree of overlap, the tolerance of such overlap will be measured in our future studies.

During precise registration, the B-spline grid size was restricted to a maximum of $22\times 22$ control points. While this cannot satisfy the accuracy requirement down to a single pixel, a larger grid size would have exponentially increased the computation cost. To account for this less-than-ideal grid size, we chose to repeatedly reuse a limited grid of control points combined with interpolation of regions showing high registration deviation, which forewent the requirement for a high sampling grid of control points. Additionally, we equally segmented the vasculature subtraction map into $10\times 10$ blocks, with each block comprising $50\times 50\text{pixels}$. The relatively large pixel size was used to make sure the selected regions would have enough information for accurate registration in the next step of the protocol.

While the protocol and accompanying data presented here are encouraging, a potential limitation for the current dual-volume scanning procedure could be that the motion-derived artifacts cannot be adequately handled when the same motions repeatedly occur across volumes. As discussed by Hendargo et al.,²³ however, this may be more of a problematic issue in ophthalmology OCT scans because subjects tend to temporarily shift their eyes when the beam passes near the periphery of the fovea, causing motion-derived artifacts to be reproducible in successive volumes. In skin OCTA, the distribution of artifact lines is random with no conceivable association to the beam scanning pattern. This possible limitation decreases with, and may be overcome by simply increasing the number of acquired volumes ($n\ge 3$). Such a method would admittedly require increased computational power, but with further refinement, this increased demand could be limited.

Additionally, a number of the methods outlined here may need reviewing if certain trends in OCT development continue. For example, the computation of current motion compensation is performed in 2-D, based on en face view of vasculatures, and the artifact lines used for frame group division are recognized from the projected OMAG signals of corresponding motion involved frames. While these are currently favorable methods, the skyrocketing speed of OCT A-scan acquisition means the OCTA algorithms are likely to rapidly change from interframe to direct intervolume.³⁹ In that scenario, the artifacts would not appear as bright lines, and accordingly, proper artifact detection methods and image registration methods would need be developed from 2-D to 3-D domains.

5.

Conclusion

We have proposed and demonstrated an automatic motion compensation method for robust motion-artifact removal in human skin OCT angiography. The compensation was achieved by sequentially performing rigid and nonrigid image registrations on dual-volume datasets to form an artifact-free and quality-enhanced vasculature composite. In general registration, a rigid affine registration was first performed on large vessels in the deep reticular dermis. The affine registration parameters were then applied to the all-depth vasculatures with a subpixel registration being added to refine the process ascribed to the involvement of smaller superficial plexus. Then, a precise registration algorithm based upon B-spline free-form deformation was used to adjust regional mismatches, which was feedback controlled through a loop of vasculature evaluation and regional interpolation. The final vasculature composites have shown considerable signal-to-noise improvement with substantial motion-artifact removed. We expect that this method would have great potential to improve the clinical diagnosis of OCT angiography in dermatology applications.