Time-space Fourier κω' filter for motion artifacts compensation during transcranial fluorescence brain imaging

Over the last decades a considerable attention has been given to optical-based imaging microscopy modalities especially to the intravital microscopy (Hillman et al 2011). To enhance the quality of microscopy images in visible and near infrared spectral range contrast materials are very often in use (Kalchenko et al 2010, Kalchenko et al 2011). Optical imaging has an ability to acquire data at high speeds; a feature that enables it to not only observe static distributions of contrast, but to probe and characterize dynamic events related to physiology (Kalchenko et al 2012, 2014, 2015, 2019a, 2019b), disease progression and acute interventions in real time (Hillman et al 2011). For example in vivo intravital microscopy (IVM) has become the mainstream technology for life science and it is used in various ways in order to extract quantitative information about essential parameters including cell location, cell motility, cell interactions or blood flow (Kalchenko et al, 2010, Hillman et al 2011, Kalchenko et al 2011, 2012, 2014, 2015, 2019a, 2019b).

One of the challenges in all live organisms imaging studies is related to the compensation of motion artifacts. Those motion artifacts should be compensated using active control approach (e.g. mechanical stabilization or gated/triggered image acquisition schemes) (Vinegoni et al 2015) and if it not possible should be removed using various types of computation (Dunn et al 2014). Most of the available computation solutions are dedicated to solve this problem and can do it with different levels of success (Dunn et al 2014, Vinegoni et al 2015). However, the changes in the illumination distribution caused by the fast dynamic changes of the forms and patterns over time due to the different flow rate of the fluorophores in veins and arteries may occur. In this way, the available methods are not efficient in cases where motion artifacts are accompanied with strong intensity distribution changes inside the region of interest during image acquisition (Kalchenko et al 2012, 2014, 2015, 2019a). It is especially relevant to the recently introduced technique (by our group), the Transcranial Optical Imaging (TOVI). TOVI is dedicated to use dynamic fluorescence as one of the contrast for imaging characterization and quantification of cerebral blood vessels through the intact cranium (Kalchenko et al 2012, 2014, 2015, 2019a). This fluorescence imaging technique is also characterized by fast dynamic changes of the forms and patterns over time related to the flow of the fluorescent material in the vessels. Nevertheless, natural movements made during the image recording produce distortions that are unique in each frame. These motions include rapid jerks or saccades, slower drifts, and high frequency tremors perturbing the temporal analysis.

Image registration is the process in which two or more images acquired at different times, from different sensors with different resolutions or dimensions, or from different viewpoints/perspectives are matched to one another (Brown 1992, Dunn et al 2014). This matching is accomplished via a process in which all of the images in the data set are transformed and aligned into a shared coordinate system. It is often described as finding an explicit function that performs a mapping of a target image onto a source image (Zitova and Flusser 2003). Also, variations in the image has been classified by three major types for the better determination of appropriate transformation (Brown 1992). The first type is the variations due to the misalignment removed by a spatial transformation when registering. Since the application of an optimal transformation in this class will remove these distortions, the variations are called distortions. The second type of variations are also distortions but are not corrected by the registration transformation because that affects the intensity values. They may also be spatial, such as perspective distortions. Finally, the third type is called variations of interest, corresponding to the differences between the images that may be spatial or volumetric but are not removed by registration. The third type variations are due to scene changes for which it is not possible to use usual algorithms, because no easy-to-find landmark locations in the image can be found along the time sequence.

Our motivation is to introduce a relevant method that could be applicable for the stabilization of intra-vital images characterized by fast dynamic changes of the illumination forms and patterns and significant motion artifacts during acquisition with variations of second and third types. The uncorrected distortions as well as the variations of interest, which together are called uncorrected variations, are not removed by registration since an exact match is not possible. The distinction between uncorrected distortions and variations of interest is important, especially in the case where both the distortions and the variations of interest are local, because the registration method must address the problem of removing as many of the distortions as possible while leaving the variations of interest intact (Zitova and Flusser 2003).

The Fourier shift theorem was proposed for registration of translated images (Bracewell and Bracewell 1986). It computes the cross-power spectrum of the sensed and reference images and finds the location of the peak in its inverse. Fast Fourier Transform can be used for shift detection in 2D and 3D spaces. It has been demonstrated that the solution to this constrained optimization problem is given by the Lagrange multiplier method (Araiza et al 2006). The method shows robustness when time varying illumination disturbances occur (Stone et al 2001). The Fourier Mellin transforms are also used to register the images, which are not only translated, but also rotated and with a change of scale (Božek and Pivarčiová 2013).

Herein, we proposed an adaptive time-space Fourier method (the κω' Fourier filtering) that was initially designed to remove the effect of intensity oscillations of astronomical observations (Tarbell et al 1989) and, to correct image motions at the limit of detection due to the noise. We also demonstrate that this approach is able to correct the entire image sequences without supervision and assignment of a reference image. In section 2, we present briefly the experimental setup to collect a fluorescent temporal sequence. The κω' method is described in section 3. Current Image Quality Assessments (IQAs) are difficult to interpret in the particular context of fast dynamic intensity associated to versatile appearance of structures in the region of interest. Quantitative assessments using both artificially corrupted images and real microscopy (Ferrari et al 2015) or using a dissimilarity index calculated between a reference image and each image in the sequence in order to detect the presence of artifacts have been recently proposed in intra-vital video microscopy (Soulet et al 2013). We introduce a useful temporal code color-coding method to perform a quick evaluation of the image quality of the entire sequence in section 4. The obtained fluorescent images collected during transcranial mouse brain imaging are shown in section 5, whereas the discussion raised by the results of this study and other practical implications of the approach are represented in section 6.

A standard fluorescent zoom microscope SZX12 RFL2 (Olympus, Japan) coupled with a CCD camera Pixelfly QE (PCO, Germany) are used for the image acquisition. Camera control and image acquisition are performed through CamWare software (PCO, Germany). A standard fluorescent illumination source was used, namely, a mercury short-arc discharge lamp. The local IACUC committee approved all experimental procedures. The Institutional Animal Care and Use Committee (IACUC) approved all animal procedures. Anesthesia was performed by intraperitoneal injection of ketamine (10 mg kg⁻¹) and xylazine (100 mg kg⁻¹) mixture, as described previously (Kalchenko et al 2014). BALB/C mice were used in this pilot study in order to verify the usefulness of the proposed method. As it was also previously described that following administration of general anesthetics, an initial cut was made, and the skin over the frontal, temporal, occipital and parietal regions was removed by blunt dissection (Kalchenko et al 2014). The exposed area was constantly moistened with saline. The mouse was then placed under the microscope lens on a special mouse holder with a warming plate, which maintained constant body temperature of 37 °C and other vital signs. For fluorescence imaging, a dose of 0.01 mg fluorescein in a volume of 50 μl was injected into the tail vein. Experiment duration was less than 1 h, after which all animals were sacrificed by barbiturate overdose. After contrast material administration, 800 raw images (exposure time: 45 ms per frame) were acquired.

The main objective of the proposed research is to remove as many of the distortions as possible while leaving the variations of interest intact in order to envision a global mapping of the brain vascularization. For instance, the temporal distribution of a fluorescent material through the vasculature is different for veins and arteries. Figure. 1(a) demonstrates color-coded image representing frame of the maximum fluorescence in each image sequence (first 300 images from the acquired stack) pixel during mouse brain visualization. The reddish colors correspond to the arteries with the fastest contrast agent administration. The greenish colors correspond to the veins which has slower contrast agent administration ability comparing to the arteries. The bluish colors correspond to the vessels with the slowest contrast agent administration (e.g. sagittal sinus). Thus, it is clearly seen that different vessel types are not observed at the same time and are localized in different places, so that a reference image to apply usual registration do not match the entire sequence when trying to compensate the disturbing motions. There exist computational imaging algorithms that are able to register sets of images, even for super-resolution, but are inaccurate when dealing with fast dynamic changes of the illumination forms and patterns and significant motion artifacts during the fast dynamic data collection (Elad and Hel-Or 2001, Nguyen et al 2001, Catlin and Dainty 2002, Farsiu et al 2004).

Zoom In Zoom Out Reset image size

We propose to apply global motion compensation on the entire sequence using a time-space Fourier method, initially designed to remove the effect of intensity oscillations of astronomical observations (Tarbell et al 1989, Molodij et al 2010). The method consists in the determination of a three dimensional filter to adjust on the Fourier spectra of the raw fluorescent image sequence. Let us consider a single three-dimensional function of intensity I_x, I_y for space, and time t that defines the fluorescent image sequences, with the notation r = (x,y) and κ = (κx, κy). The filtered time-space intensity is applied in the κω' space where κ and $\omega '$ are the spatial and temporal frequencies. Let $\tilde I = \left( {\boldsymbol{\kappa},\omega '} \right)$ be the Fourier transform for space and time of the Intensity I(r,t). In the spatio-temporal analysis assumption, the Fourier transform of the fluorescent signal I(r,t) consists of a projection of the temporal signal onto the orthogonal basis $ex{p^{ - 2i\pi \left[ {\omega t} \right]}}$ and the spatial signal onto the orthogonal basis $ex{p^{ - 2i\pi \left[ {\mathbf{r}\boldsymbol{\kappa} } \right]}}$, presented as:

$$\begin{equation} I\left( {\boldsymbol{\it{r}},t} \right) = {\rm{ }}\mathop \int\!\!\!\int \limits_{{R^2}}^{} \mathop \int \limits_R^{} {\tilde{\!I}} \left( {\boldsymbol{\kappa},\omega '} \right)ex{p^{ - 2i\pi \left[ {\boldsymbol{r\kappa} + \boldsymbol{\omega} t} \right]}}W\left( {\boldsymbol{\kappa},\omega '} \right)d{\boldsymbol{\kappa}}d\omega'{\rm{ }} \end{equation} \tag{ 1 }$$

where the 3D time-space window filter $W\left( {\boldsymbol{\kappa},\omega '} \right){\text{ }}$ is defined as

$$\begin{equation}W(\kappa ,\omega ') = \left\{ {\begin{array}{*{20}{c}} {1,{\rm{ }}if{\rm{ }}\sqrt {\kappa _x^2 + \kappa _{y{\rm{ }}}^2} \frac{{\omega '}}{V}}\\ {{\rm{0, }}elsewhere} \end{array}.} \right.\end{equation} \tag{ 2 }$$

Figure. 1(b) shows an illustration of the reciprocal time-space representation of the Fourier analysis of the fluorescent signal. Fourier time-space is represented taking into account the periodicity and symmetry properties of the Fourier transform; $\omega '$_max appears twice while $\omega '$_min is indicated at the center. The total sequence duration and the time sampling determine, respectively, $\omega '$_max and $\omega '$_min. The pink cones delimit the removed frequencies of the Fourier space-time, defined by the velocity $V = \omega '/k$. In the presented example v = 3.10^–3 ms⁻¹.

The use of a local spatiotemporal analysis allows the discrepancy of different regions in the Fourier space. For instance, the time-space filter of equation (2) defines a cone of velocity $V \approx \omega '/k$ that determines the velocity of the motions to discard. All Fourier components inside the κω' volume are removed. Note that this velocity parameter V is a mathematical parameter of the underlying geometric transformation law between space and time. The attempt to relate this velocity to a physical process requires a more difficult and dedicated experimentation procedure.

In figure 2, Filters are displayed for different values of the ω'/κ ratio. The method acts as a selection of frequencies distributed in the 3D volume of the Fourier Space.

Zoom In Zoom Out Reset image size

In figure 2 a linear relationship between space and time has been assumed to determine the values of the frequencies to remove in the κω' space. Another quadratic behavior might also be envisaged in the case of a diffusive law assumption to take into account the localization of the vessels at different depths.

As it was mentioned before, it is not possible to use a reference image for motion artifacts compensation as far as fast intensity distribution changes such as arteries and veins are collected at different times and located at different places. Figure 3 shows the result of the processing after applying a stationary signal detection utilizing a reference image by comparison to the κω' method. The rigid registration method computes the cross-power spectrum of the sensed and reference images and finds the location of the peak in its inverse (implemented in the furnished Matlab code for comparison). The rigid registration method produces saccades due to the dissimilarities between the selected reference and the current image as indicated in figure 3(a) where is displayed the blurring coming from the misalignment after the registration.

Zoom In Zoom Out Reset image size

As a practical tool to obtain a quick look of the efficiency of the described method, we propose to use a temporal code color to overcome the difficulty due to the blending of the motion effects and the change of structures along the time. The main idea is to weight each image with a different color before the average in order to regroup structures appearing with the purpose to distinguish the change of structures from the motions. Color-coded maps are created using a multicolor merge approach of a full-unprocessed time sequences using the temporal color code plug-in in Fiji image analysis software (http://fiji.sc) (Schindelin et al 2012). Briefly, 300 frames are color coded and merged for a qualitative representation of image/time distortions. A 'rainbow RGB' color scheme is chosen to correspond to the different levels of intensity distribution appearance along the temporal sequence.

By comparison, we evaluated the image quality using methods such as the Similarity Image Assessment (SIA) (Molodij et al 2014), the entropy (Kersten 1987, Molodij et al 2015), and the sharpness (Molodij et al 2015) on the final obtained average image to compare with the time color image distortion method. This effort was first directed toward tuning quickly the $\omega '/k$ parameter to find the optimum parameter. The lack of a common reference that fits the entire sequence leads us to adopt absolute assessments such as the entropy to estimate the quality of information retrieved or from the sharpness metric using the knowledge provided by adaptive optics in astronomy (Muller and Buffington 1974, Zommer et al 2006, Molodij et al 2015). Sharpness assessment characterizes the features in the image while the entropy method indicates the level of information in each frame. An advantage for the entropy assessment is coming from its relative insensitivity to the intrinsic changes of the features. Optimum value of the filtering velocity parameter V can be derived and eventually automated.

Moreover, we remarked that the κω' filter acts to remove the 'background noise' due to the fluorescence scattering in the medium for large value of the velocity. As an example in figures 4(b) and (d), hidden structures emerge for recognition when applying the proposed filtering. Nevertheless, the intensity of the pixels are modified and maybe not realistic. An accurate study of the image quality assessments to quantify the reliability has to be studied carefully and is beyond the scope of this paper. In regard only of the IQAs entropy scores and denoting that the κω' method depends on the temporal size of the sequence, we propose to adapt the κω' depending on the sequence time duration.

Zoom In Zoom Out Reset image size

As one can see in figures 4(a) and (c) the κω' filtering is not applied when $\omega '$=${\text{ }}\omega '$_min in order to keep the 'background noise'. In the case of a longer temporal sequence (as indicated in figure 5), the processing leads to a similar result when determining the best filter, i.e. one obtains the same value of the optimum velocity from the IQAs.

Zoom In Zoom Out Reset image size

In the case of short temporal sequences, applying a rigid registration eventually followed by the κω' would be more adapted, when dynamic changes of the illumination are less prominent. Figure 5 plots the IQAs values for a long temporal sequence case study (the time sampling is longer than the spatial dimensions). For the sake of simplicity, the maximum score to display the values on the range from 0 to 1 normalizes all values. IQA indicates a value of V = 0.9 as the optimum value for both types of filtering. The sharpness of the long duration sequences is significantly decrease at the value of V = 3. Important to notice that 3D Fourier transform computed with indirect dimension (2 spatial and 1 temporal). Depending on the size of the 2D spatial image (detector size), the time sequence has to be at least the same dimension. Long temporal sequences are defined by a time resolution larger than the size of spatial x or y-axis. A value of the velocity of one of the κω' method is often the best solution for a cube of size x = y = T.

We verified the validity of the linear relationship assumption for the κω' filter by increasing the duration of the sequence (from 300 up to 800 frames). The corresponding optimum value of the velocity V is then proportionally shifted.

In the experiment related to figure 3, the tissue was grossly immobilized, leaving only minor distortions in the field resulting from respiration and heartbeat. However, in some cases, a motion can actually shift the entire field of view. In this close-up image, the temporal sequence was made up 800 images and the growing factor scale was about 20 by comparison to the following experiment. There is an interest for an accurate compensation of jerks before any characterization or quantification of cerebral blood vessels through the entire intact cranium shown in figure 6.

Zoom In Zoom Out Reset image size

This experiment was also an opportunity to apply the κω' method on an irregular data cube (x = 696 pix, y = 512 pix, t = 300 × 45 ms). Obviously, a better spatial resolution is obtained when applying the filtering. But important also, is the accurate recognition of the structures that appear at different moment and are grossly indicated by the color-coding. For instance, in figure 6(a), due to the involuntary motions, the discrepancy between veins structures (indicated in blue and green) at different instant becomes difficult. By comparison, in figure 6(b), one can expect to label the different functional vessels. In the κω' Fourier space mixed information is coming from image structures as well as dynamics. The velocity parameter V delimits the Fourier time-space corresponding to the undesirable motions faster than the time dynamic scale of the illumination changes. If the temporal sampling rate of the sequence is high enough, tuning the velocity parameter can separate both effects. Saccades from jerks are described by velocities larger than characteristic velocity of illumination pattern modification in the image, are easy to extract.

Figure 8 displays the time color image distortion of the fluorescent sequence when increasing the value of the velocity. Image quality entropy assessment indicates a value of V = 0.9 as the optimum value for the registration (see figure 7). Important to notice that sharpness of the long duration sequences is significantly decrease for the value of V > 3. Nevertheless, increase the number of frequencies beyond the optimum to remove totally the jerks is accompanied by a loss of information in relation with the image patterns.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The high spatial resolution structures vanish and complex structures are mixed. Smaller looking vessel diameter for larger velocity parameter is due to some 'erosion effect' related to the loose of corresponding frequencies. Images look sharp but we remark that details are vanishing. This effect may explain the relatively low scores of both entropy and sharpness IQAs for large value of the velocity in figure 7. It is important to notice that both sharpness and entropy IQAs are applied directly on the average sequence image after κω' filtering without any correction of the dynamic of the signal that is incompatible with the definition of the entropy IQA (see the Matlab code, available at {stacks.iop.org/PMB/65/075007/mmedia}).

In this section, we discuss first issues related to the definition of the registration. The novelty here is that we apply the algorithm in a 3D space. Many of the problems arising from the projection of 3D space onto a 2D image may no be longer relevant (Brown 1992). From 2D studies, it has been accepted that the fundamental characteristic of any image registration technique is the type of spatial transformation or mapping used to properly overlay the images (Brown 1992). In this paper, our motivation was to introduce a relevant method to remove only some of the variations; the effects of illumination changes may be difficult to remove, or we are not interested in removing them, i.e. there may be changes that we would like to detect. So, the distinction between uncorrected distortions and variations of interest is important, especially in the case where both the distortions and the variations of interest are local. Because the registration method must address the problem of removing as many of the distortions as possible while leaving the variations of interest intact, it becomes necessary to distinguish between whether certain variations are global or local and whether the selected transformation is global or local. For example, images may have local variations, but a registration method may use a global transformation to align them. In the κω' method, the difficulty to display an explicit geometric transformation comes from the difficulty to distinguish the frequencies related to the global and local variations in the κω' space. Only an underlying geometric transformation law can be suggested.

In the purpose to determine the limitations and optimize the parameters, we applied the κω' processing on synthetic data. In the following simulation the results of which are shown in figure 9, the field of view is parsed with three-colored number of different sizes.

Zoom In Zoom Out Reset image size

Each color corresponds to a specific variation of the intensity indicated in figure 9(b). We simulated distortions such as respiration and heartbeat together with random jerks. Plots of the generated motions are indicated in figure 9(c). The noise and the scattering effect are not generated. This simulation addresses the following purposes, in order to:

• Define spatial patterns for recognition (i.e. the three numbers) that are associated to three different time behavior (such as the brain vascularization);
• Produce saccades when applying a rigid registration (generating enough dissimilarities between the reference and the current image);
• Mix frequencies components close to the real case (such as respiration and involuntary random jerks);
• Blend the information coming from the intensity variations and displacement motions in the Fourier space (that is not the case for a simple moving disk, for instance);
• Show the interest of the time color image distortion method.

In figure 10 simulated images before and after the κω' method are displayed using the average time color image distortion. The entropy IQA indicates a value of V = 0.4 as the optimum value for the correction in agreement with the subjective quality from the time color. The image looks sharper but smallest details are not resolved. In the following, we show that this residual blurring is coming principally from the large slow motions characterized by low velocities close to the intensity light variations, fastest jerks being corrected by the κω' method.

Zoom In Zoom Out Reset image size

In the following, we simulate the jerks as a simple sinusoid defined by A*sin[2πF*t] where A is the amplitude and F the frequency of the sinusoid. We keep the three different time variations of intensity indicated by both the three color and numbers. Figure 11(a) shows the reference sequence image sharpness case when no motions are simulated is displayed. Figure 11(b) shows the subfield indicated by a red square for the purpose of a close-up comparison.

Zoom In Zoom Out Reset image size

Figure 12 presents a comparison of the time color image distortion before (a), (c) and after the κω' filter (b), (d) for two simulated sinusoidal motions A*sin[2πF*t]. It was found that the κω' method is first affected by the frequencies and at a second order by the amplitude. For instance, the processing performs identically for motion amplitudes smaller than 15 pixels when the frequency of the oscillation is 1/300. Entropy IQA indicates the same optimum value for the κω' filter velocity (a value of V = 2). By comparison, in the case of a smaller frequency with a smaller amplitude, the more difficult is to compensate the motion closer such as the oscillation with A = 5 and F = 1/50. The discrepancy in the Fourier space between motions and the variations of intensity is more difficult to determine.

Zoom In Zoom Out Reset image size

Figure 13 shows the corresponding subfield close-up for the two simulated cases. In this simulation, the sinusoidal frequencies refer to the collection duration time. Longer is the temporal sequence; more efficient is the κω' method.

Zoom In Zoom Out Reset image size

Despite the fact that it is not possible to display a geometric transformation map, the proposed method addresses the problem of removing as many of the distortions as possible while leaving the variations of interest intact. At least, usual registration methods face difficulties to correct occasional images motions when jerks are mixed to the dynamical changes. Image registration technique such as the local image correlation leads to an unsatisfactory result due to the difficulty to determine a static reference all along the temporal sequence. Attempts to determine the relative motion of each image in respect with the previous one in the temporal sequence leads to introduce large saccades that compromise the patterns recognition.

Another problem which can be potentially solved by the proposed method is removing of artefacts caused by the highly scattered light from deep layers of tissue. The possible way of image quality improvement in this case can be application of a diffusion law instead of a velocity in the reciprocal Fourier space for instance.

In the case of short temporal sequences or at the end of a fluorescent sequence when the dynamic slows considerably, we applied successfully a rigid registration followed by the κω' filter to improve the spatial resolution. The κω' method is able to preserve the movements of tiny illumination patterns of interest, and complete the correction of usual linear registrations.

We presented a time-space κω' Fourier approach for the stabilization of intravital fluorescence-based imaging in the presence of fast dynamic changes of the observed imaging patterns, and significant motion artifacts during their acquisition. The introduced technique is an alternative to correction of variations in intravital images of brain obtained through the intact cranium in vivo. Despite fast changes of the intensity patterns with time, the current approach is able to improve the spatial resolution of the entire sequence, and eventually to reduce the 'background noise' emerged during the fluorescence imaging. The proposed adaptive approach would be highly beneficial when combined with standard rigid registration or deconvolution methods.

We also proposed a simple procedure, based on the time visual image distortion method, for a quick estimation of the compensation quality. The proposed technique can be used autonomously without supervision and assignation of a reference image, the ideal velocity threshold being determined by testing different velocities and comparing the entropy.

We developed codes for both the Matlab and IDL (Interactive Data Language) environments with a similar result although the Fast Fourier Transforms procedures are slightly different. We furnished the Matlab codes of the adapted κω' method depending on the duration of the sequence or/and including a rigid registration step. The time cost to process a sequence of 600 frames is less than a minute for a current laptop working with Matlab.

The authors would like to thank The Henry Chanoch Krenter Institute for Biomedical imaging and Genomics for the main support of this work (Dr. V Kalchenko, 'Weizmann Staff Scientists Program'), and support provided by COST CA16118—European Network on Brain Malformations. This work has been also supported by the European Union's Horizon 2020 research and innovation programme under grant agreement No 863214. Professor Meglinski also acknowledge partial support from Academy of Finland (project 326204), MEPhI Academic Excellence Project (Contract No. 02.a03.21.0005) and the National Research Tomsk State University Academic D.I. Mendeleev Fund Program. Anton Sdobnov acknowledges the Finnish Cultural Foundation (00180998) grant. We also thank reviewers for valuable comments.

The presenting authors do not have any financial conflicts of interest regarding the content of this manuscript.