The results of this study show that our algorithm was effective in suppressing decorrelation caused by cardiac and respiratory motions of the animal. Cardiac motions mainly caused GPFs, while respiratory motions caused BISs. Subsequent application of the GPF, BIS, GPF and NGPF corrections, the best combination among the various possibilities, significantly reduced the decorrelation, thus maintaining a high average image correlation of 0.76 ± 0.04 even after three seconds measured from five animals (). In addition, the combination remarkably stabilized the phase of the OCT signal ().
4.1 Choice of the reference frame
The choice of reference frame affected the performance of motion correction processing. This study chose the reference frame where cardiac and respiratory motions were most minimized. When we chose a reference frame from the middle of the cardiac/respiratory motions, the cross-correlation of raw data was lower and noisier. The mean cross-correlation magnitude of raw data decreased to 0.25 from 0.60. Also, the performance of motion correction was worse than that presented in this paper. The mean cross-correlation magnitude of the motion-corrected data was lower, 0.76 as opposed to 0.79.
Although we chose the first frame as the reference to better show decorrelation in time, it would be better during practical processing to choose a frame from the middle of the total acquisition time. This would lead to smaller noise because the performance of BIS and GPF correction is degraded as the time gap between the two frames increases. We used an algorithm to automatically find an optimized reference frame: every frame during one cycle of the cardiac motion in the middle of the acquisition time was tested to find the one producing the highest mean cross-correlation of raw data.
4.2 The number of upsampling
We upsampled images in the BIS correction to find sub-pixel shifts. This upsampling, however, requires a large amount of computation, thus making the BIS correction take much longer time than the GPF correction. Further, the computational load increases by the fourth power of the amount of upsampling. We tested various upsampling numbers from 2 to 8 and chose 4 because GPFs, which are partially coupled to sub-pixel BISs, were additionally corrected.
As an alternative means of compensating for sub-pixel shifts without the large computational load, we tested a Fourier transform-based method. First, we upsampled the cross-correlation for various pixel shifts (
Eq. (2), rather than upsampling every image. From the upsampled shift-correlation relation, we found the sub-pixel shift maximizing the cross-correlation. Then, we compensated sub-pixel shifts with the Fourier transform.
where

is the Fourier-transformed OCT signal and F
−1[ ] means the inverse Fourier transform. The computational load of this alternative method was much lower than that of the method based on explicit upsampling of every image (~1/100 smaller computation time); however, the performance was significantly lower, 0.768 ± 0.025 as opposed to 0.800 ± 0.018, in terms of the mean magnitude of cross-correlation.
4.3 Coupling between image shifts and global phase fluctuations
As described in the Results section, respiratory motions mainly caused BISs and cardiac motions caused GPFs. This relationship may not always be correct, however, as cardiac motions also can cause image shifts when they are so large as to cause supra-pixel movements of the brain. Therefore, it is more reasonable to understand that the BIS correction is more effective in compensating for motions larger than the pixel width, while the GPF correction is more effective in compensating for sub-pixel motions. Of course, the BIS correction could also suppress sub-pixel motion-oriented decorrelation if the amount of upsampling is sufficiently large. We confirmed that the amplitude of cardiac motion-oriented decorrelation () was further reduced when the BIS correction was applied with 16-fold image upsampling. This result implies that BISs and GPFs are coupled with one another, at least partially. Of course, the 16-fold upsampling required too much computational load to be of general utility.
4.4 Motion correction with a mask
As can be seen in the noise map (), noise from vessels were much larger than those from stable tissue. Furthermore, when we looked at time courses of phase signals from vessels, the noise was too irregular to be related to global fluctuations. Those irregular noises might be attributed to movements of red blood cells. These movements (1-10 mm/s) are so fast with respect to our spatiotemporal resolution (3.5-7 μm and 4 ms) that OCT signals from vessel voxels may be totally uncorrelated after 1-2 time steps. As a result of this large decorrelation the complex-valued OCT signal appeared to be moving stochastically in the complex plane. Therefore, this irregular noise from vessel voxels might contaminate the cross-correlation in both BIS and GPF corrections (
Eq. (2) and
3).
Another source of contamination was noise from the region of air and deep tissue. Since the amplitude of the OCT signal from those regions was so small, data points for the OCT signal looked like a Gaussian distribution in the complex plane. For this reason, the phase of those regions was very unstable and irregularly varying, thus contaminating the cross-correlation.
In order to minimize these contaminations, we examined a method using a mask. First, we applied the C3 method to find a relatively accurate noise map, and built a mask based on both noise and intensity maps. The mask had a value ranging from 0 to 1 at each voxel depending on its noise and intensity. Then, we used the mask in calculating the cross-correlation (
Eq. 2 and
3) during application of the C5 method to raw data. It resulted in slightly lower (2%) noise in the tissue area (at a depth of 150-600 μm from the cortical surface). Of course, the use of a mask doubled the computation time, which can be a drawback in a practical analysis.
4.5 Feasible applications
Our motion correction processing can be helpful not only in noise reduction but also in enhancement of vessel structure visualization. Since the angiogram has been widely used in revealing vessel structures from dynamic OCT imaging data [
6], this study compared angiograms of raw data and motion-corrected data where the C5 method was used. The angiogram was obtained with the time derivative of the OCT signal.
This angiogram differs from the noise map in that it is proportional to the mean displacement of the OCT signal during one time step in the complex plane, whereas the noise map is the radius of the OCT signal distribution for the total acquisition time. For example, when a unit-magnitude OCT signal rotates π/100 at each time step for 200 time steps, the angiogram intensity is π/100 whereas the noise map intensity is 1. Therefore, the angiogram may better represent blood flow-oriented decorrelation when overall phase noise is large, because such a large phase noise results in the OCT signal rotating more than 2π for a given time (, for example).
As can be seen in
, cross-sections of vessels were much more clearly evident in the motion-corrected angiogram. In contrast, the noise map and angiogram of raw data did not show much difference from the intensity map. This is likely due to large GPFs. A large GPF causes the OCT signal to rotate more than 2π in the complex plane and, as discussed above, such rotations make the deviation of signal similar to its magnitude. The angiogram can also be contaminated by large GPFs such as these because they make the displacement of the OCT signal in the complex plane much larger than that due to blood flow-oriented decorrelation.
Although we used as an example the cross-sectional angiogram obtained from the dynamic imaging data, general angiograms do not require such long-term dynamic imaging. Further, the proposed algorithms are intended for dynamic OCT imaging, not for structural imaging including angiograms. Nevertheless, it is worth showing how the proposed method works for 3D structural imaging (angiograms, for example). We performed another OCT scan for an angiogram, where B-scans were repeated two times for each cross-sectional plane (i.e., two time-points dynamic imaging) [
6]. We used
Eq. (9) to reveal vessel structures and performed maximum intensity projection through the z axis to obtain an
en face image. As can be seen in , we confirmed that the proposed method works very well not only for dynamic imaging but also for 3D angiography.
The motion correction algorithms presented in this paper would also be helpful in studies investigating time-varying physical quantities of tissue with phase-resolved optical methods. For example, the combination C5 is currently being used in one of our ongoing studies looking into depth-resolved hemodynamic responses of the cortex associated with functional activation. The combination C5 is effectively reducing cardiac and respiratory motion-oriented noise in our functional studies, thus enabling much clearer contrast in the hemodynamic responses.