Image-based flow decomposition using empirical wavelet transform

We propose an image-based flow decomposition developed from the two-dimensional (2D) tensor empirical wavelet transform (EWT) (Gilles 2013). The idea is to decompose the instantaneous flow data, or its visualisation, adaptively according to the averaged Fourier supports for the identification of spatially localised structures. The resulting EWT modes stand for the decomposed flows, and each accounts for part of the spectrum, illustrating fluid physics with different scales superimposed in the original flow. With the proposed method, decomposition of an instantaneous 3D flow becomes feasible without resorting to its time series. Examples first focus on the interaction between a jet plume and 2D wake, where only experimental visualisations are available. The proposed method is capable of separating the jet/wake flows and their instabilities. Then the decomposition is applied to an early stage boundary layer transition, where direct numerical simulations provided a full data-set. The tested inputs are the 3D flow data and its visualisation using streamwise velocity&{\lambda}2 vortex identification criterion. With both types of inputs, EWT modes robustly extract the streamwise-elongated streaks, multiple secondary instabilities and helical vortex filaments. Results from bi-global stability analysis justify the EWT modes that represent the streak instabilities. In contrast to Proper Orthogonal Decomposition or Dynamic Modal Decomposition that extract spatial modes according to energy or frequency, EWT provides a new strategy as to decompose an instantaneous flow from its spatial scales.


Introduction
The age of big data witnesses an expeditious generation and accumulation of flow data, both numerically and experimentally (Duraisamy et al. 2019). Even standing on top of data, on most occasions, it remains not a simple task to understand a fluid flow due to the broad range that the flow may comprise in its temporal and spatial scales. On the other hand, complete data-set may not always be available. The extraction of concerned knowledge from projected flow data (e.g. from flow visualisations) represents immense potentials for practical applications (Raissi et al. 2020).
The development of model reduction and modal analysis techniques has played an essential role to promote this effort (see reviews by Taira et al. 2017Taira et al. , 2020. Figure 1 provides a general view of the well-developed data-based modal obtained by physically interpreting the wavelet coefficients. In the examination of fluid flow data, wavelet analysis has been widely used as a flexible "microscope" discriminating scales and positions simultaneously. The coherent vortex extraction (CVE) proposed by Farge et al. (2001) splits turbulent flows into coherent and incoherent parts. CVE achieves its goal by projecting the vorticity field onto an orthogonal wavelet basis, and a threshold on the resulting wavelet coefficients is then specified to identify the coherent component. In the last decade, wavelet-based numerical algorithms and turbulence modelling have been developed (Schneider & Vasilyev 2010). Their strengths lie in the competence to unambiguously identify and isolate localised, dynamically dominant flow structures which improve the efficiency of the computation. EWT provides the first mathematically rigorous method to adaptively analyse a signal (Gilles 2013). The principal idea is to build a set of adaptive filter banks based on the property of the input signal itself, such that meaningful modes are extracted, each possessing a particular section of the Fourier supports. This feature matches well with the requirement to extract modes from instantaneous flow data or its visualisation.
In §2, the formulation and algorithm of EWT for flow decompositions are introduced. §3 provides two examples as EWT being applied to 3D instantaneous flow data and its visualisation, followed by the comparison with the other methods in §4. The study is concluded in §5.

Imaged-based flow decomposition using EWT
We term the complete flow data-set or the visualisation (image or video) as a function f (x, y, z; t). Here x, y, and z represent the spatial coordinates, and the optional t denotes time. Note that flow visualisations project data onto a 2D image or video: f (x, y; t). The visualisation may be qualitative and aided with auxiliary materials (e.g. smoke, sensitive coatings), be quantitatively made with state-of-art facilities (e.g. particle image velocimetry), or from numerical simulations (e.g. using λ 2 criterion (Jeong & Hussain 1995) to identify vortex structures). We decompose the flow using two-dimensional (2D) tensor EWT (Gilles et al. 2014) as summarised in Algorithm 1.
EWT defines its scaling functionφ m and empirical waveletsψ m (see Appendix A) in a normalised Fourier axis with 2π periodicity. The set {φ 1 , {ψ m } M −1 m=1 } is a tight and orthogonal frame of L 2 (R). Consequently, the energy of the signal is conserved by the extracted EWT modes. The empirical wavelet transform of a signal f (t) is obtained by (2.1) where stands for the inner product, † is the complex conjugate, ∧ and ∨ indicate the Fourier transform and its inverse.
We show in figure 2 an example of the adaptive filter bank. The signal in Fourier spacê f (ω) is shown with a white line. In this study, the Fourier supports {ω m } are adaptively determined such that local maxima off (ω) are retained in different Fourier segments (indicated with different colours). As a result, the extracted modes strategically capture components of the input signal, each standing for a section of the spectrum. The idea of 2D tensor EWT is to perform EWT successively in two directions, with the corresponding filter banks built according to the spatial-averaged spectrum. When flow visualisation is used as an input, it is essential to have the principle direction of the flow (e.g. direction of the mean flow) coincide with either direction of decomposition, such that the averaged skeleton mode shadow mode Figure 2. A sketch in the Fourier space of the signalf (ω), the detected supports {ωm}, and the filter bank {φ1(ω),ψm(ω)}. ω ∈ [0, π] and abides by the NyquistShannon sampling theorem. The Fourier supports are determined from the spectrum of the signalf (ω) such that each EWT mode amounts to part of the spectrum. The first and the last EWT modes, W1 and WM , are termed the shadow mode and the skeleton mode, respectively.
Algorithm 1: Flow decomposition using 2D tensor empirical wavelet transform Input: f (x, y, z; t) with N x , N y , N z and N t points in corresponding indices. If t (optional) is present, the input is a 3D unsteady flow or its visualisation (video); M x , M y , the number of desired filters/modes in the direction of x and y. Note that we show the decomposition algorithm in the x and y directions, while in practical applications, the two directions can be chosen from x, y and z. 1 For inputs with flow visualisations, detect the principle direction of f (x, y, z; t)(details can be found in §3.1). Rotate f such that the principle direction is in line with x or y; 2 In the x direction, compute its FFT to obtain the x-averaged spectrum. For the input with N t snapshots, the spectrum is also averaged in time: Nz k=1 Nt t=1 f (j, k, ω x ; t) ; 3 Similarly, obtain the y-averaged spectrum: Nz k=1 Nt t=1 f (i, k, ω y ; t) ; 4 Perform the Fourier boundary detection onf x (ω x ) andf y (ω y ) respectively; 5 Build filter banks B My−1 m=1 } using (A 1) and (A 2); 6 Filter f along the x direction with B x , resulting in M x outputs; 7 Filter each output in step 6 along y direction using B y ; 8 If necessary, reconstruction of the input f (x, y, z; t) is performed by inverting the y and x filtering successively following Output: EWT modes: W i,j (x, y, z; t) with i = 1, ..., M x , j = 1, ..., M y .  spectrum maximises its physical relevance. This corresponds to step 1 of algorithm 1 and will be detailed with an example in §3.1.

Applications and results
We show two applications of flow decomposition using EWT: the interaction between a 2D wake and a jet plume, and the early-stage boundary layer transition subject to free-stream turbulence. The two examples serve to demonstrate the applicability both to flow data and its visualisation, irrespective of their experimental or numerical origin.

Interaction between a jet plume and a 2D wake (images from experiments)
In this example, the decomposition is applied to the experimental observation of the interaction between a jet plume and a 2D wake (Roquemore et al. 1988). The slot jet is located in the centre of the rectangular face of a bluff body. The flow is visualised with reactive Mie scattering technique, where streamlines are highlighted by TiO 2 particles that are spontaneously formed by the reaction of TiCl 4 in the jet with water in the wake. The flow visualisation is shown in the first column of figure 3 with jet velocity increasing from 18.5 cm/s in panel (a) to 55.5 cm/s in (c). As introduced in Roquemore et al. (1988), the flow is initially dominated by the wake flow when the shear layer velocity of the jet U jet is smaller. Along with the increase of U jet in panel (b), wake instabilities are considerably prohibited, and only wavy structures are observed. In panel (c), the flow is dominated by the jet instability as indicated by the change in the direction of rotation of the vortices.
The flow decomposition is performed with M x = M y = 3, leading to nine EWT modes, W i,j (i, j = 1, 2, 3). The lower modes (columns 2 ∼ 5 of figure 3) well isolate key components of the flows. In particular, mode W 1,1 (the shadow mode), comprising the lower-end spectrum both in the flow and traverse directions, can be used to detect the principal direction of the flow. Mode W 2,1 possesses mid-spectrum in the traverse direction and features the mean flow. Representing the mid-spectrum of the flow direction, mode W 1,2 and W 2,2 present the flow instabilities with W 2,2 accounting for higher wavenumbers in the traverse direction. Comparing the three cases adheres to the intention of the experiment: an increase of the jet velocity leads to the suppression of wake instabilities (panel b), and promotion of the jet instability downstream (panel c), as it is supported by W 2,2 in the three panels. It is also helpful to inspect the flow by combining certain modes: W 1,1 ⊕ W 2,1 and W 1,2 ⊕ W 2,2 . The combination is obtained by performing an inverse EWT with corresponding modes. As it is shown, W 1,1 ⊕ W 2,1 and W 1,2 ⊕ W 2,2 account for a wider spectrum in the traverse direction and differentiate in the flow direction, thus delivering a more comprehensive view of the mean flow and wake/jet instabilities respectively. Note that we have reconstructed the higher modes (the last column of figure 3) with W 1,3 ⊕ W 3,1 ⊕ W 2,3 ⊕ W 3,2 ⊕ W 3,3 . The higher modes amount for the rest of the spectrum and retain all the small scales of the visualisation. Note that by increasing the number of EWT modes, the higher modes (the skeleton) will contain less information of the input.
Recalling algorithm 1, in 2D tensor EWT, the wavelet filter bank is built based on the averaged spectrum in the x and y directions of the signal. It is crucial that the flow's principal direction (for example, the direction of the mean flow, orientations of the geometry) coincides with these directions, such that the flow physics can be correctly separated. In the 37.0cm/s case of the above visualisation, we have rotated the source counterclockwise by 2.0 • before the EWT is applied (step 1 of Algorithm 1). We take this case as an example to show how to accurately detect the principal direction. As shown in figure 4, the shadow mode W 1,1 , which is free from small scales, characterises the principle direction of the signal. This direction is thus numerically recovered by the gradient, |∇(W 1,1 )|. The principle direction is the direction along which the average gradient reaches a minimum, as shown in figure 4(c). Consequently, it is found that the source shall be rotated counterclockwise by 2.0 • . Figure 4(d) shows the rotated image as compared with figure 4(a). Panel (e) shows the decomposition when the source image has an inappropriate orientation at which W 2,1 and W 2,2 fail to capture the mean flow and instabilities.

Early stage boundary layer transition (from direct numerical simulations)
Laminar-turbulent transition prompts a significant broadening of flow scales that ultimately takes the form of coherent structures. This example examines an incompressible boundary layer flow over a flat plate with an elliptic leading edge (see figure 5). To promote flow transition, free-stream turbulence (FST) from a nonlinear optimisation procedure (with the target that the maximum perturbation energy is reached at a designated time) is specified ahead of the leading edge. The width of the domain is chosen such that three discernible streamwise-elongated streaks are generated. A periodic boundary condition is employed in the spanwise direction. Further numerical details and set-up of the simulations can also be found in Wang et al. (2019). In this example, we first apply the algorithm to flow visualisation, followed by the decomposition with instantaneous 3D velocity data. The flow has been visualised through its pivotal structures, i.e., with iso-surfaces of λ 2 = −0.02, u = 0.8 and contour slices of u = 0.1, 0.2, ..., 0.9 as presented in figure 5, here u stands for the streamwise velocity (normalised by the freestream value). A grey-scale top-view image with the same iso-surfaces serves as the input of the image-based flow decomposition. As can be seen, the 2D image features the generation of streaks and the formation of hairpin vortices.
We show the flow decomposition in Figure 6(a). In this example, we have M x = M z = 3. By construction, mode W 1,1 (shadow mode) and W 3,3 (skeleton mode) represent the flow structure with a spectrum that amounts to the lower-and higher-end both in the streamwise and spanwise directions. Mode W 1,2 and W 1,3 isolate the streaks. The secondary instabilities of streaks are captured by modes W 2,1 , W 2,2 and W 2,3 , featuring the low-to-high wave numbers in the spanwise direction. From W 2,1 , it is seen that Streak A gives rise to sinuous instabilities after x = 40, followed by the instability of Streak C. It also indicates that perturbations around Streak C have a larger amplitude after x = 70. These observations are however hidden in the source image. Finally, in mode W 3,1 and W 3,2 , localised helical vortex filaments are identified, marking the meandering of streaks and appearance of smaller scales.
To verify the observation from EWT modes, a biGlobal stability analysis (secondary instability of streaks) is performed at x = 43, 44, 45 (to secure the onset location of secondary instabilities) and x = 50, 70 (to identify multiple modes). We plot temporal .., 0.9 in cross sections that are perpendicular to the freestream direction. Both the iso-surfaces and contours are coloured by u. We apply imaged-based flow decomposition on the grayscale 2D image that contains the same iso-surfaces. The image has a top view on which the flow data in the wall normal direction is projected. Within the computation domain, three discernible streaks are generated, and they become unstable downstream before giving rise to hairpin vortices.
growth rates versus streamwise wavenumber α in panel (b). Their peak growth rates are marked with red crosses, and the corresponding eigenfunctions & base flows are provided in panel (c). A positive growth rate is seen at x = 44, which demonstrates the onset of secondary instabilities. This mode (termed mode 1) has its eigenfunctions localised around Streak A. Further downstream, the flow is unstable to multiple modes.
From around x = 50, mode 2 (centred around Streak C) also becomes unstable, though, at a smaller growth rate. However, further downstream at x = 70, mode 3 (centred around Streak C) becomes the most energetic. All these secondary instability modes are of sinuous nature. Generally, the stability analysis matches rather well with EWT modes; the onset of multiple secondary instability modes, their sinuous nature and the dominance of instabilities around Streak C are evidenced by mode W 2,1 . For the case of a video input, an unsteady flow decomposition is provided as a supplementary file.
The video records the flow from t = 0 when FST is introduced to the laminar flow until the time step illustrated in Figure 5 and 6. The video-based EWT modes provide the evolution of flow structures subject to the temporal-spatial averaged Fourier supports.
To test the sensitivity and robustness of the method with respect to the input image, a comparison of EWT modes with different flow visualisations is shown in figure 7. Images with various isosurfaces of u and λ 2 are adopted as shown in panel (a). The same filter banks from EWT are used to process these inputs. As can be seen from panel (b) and (c), the streaks and their secondary instabilities are correctly extracted by mode W 1,2 and W 2,1 for most of the cases, regardless of the different inputs. With isosurfaces of u alone, secondary instabilities are captured, but the importance of Streak C is not revealed. For a more inimical input (u = 0.5, λ 2 = −0.10), the EWT modes still result in satisfactory decompositions. We further apply EWT to the decomposition of 3D instantaneous data, as shown in figure 8. The flow has been decomposed in the x and z directions with M x = M z = 3 based on the streamwise velocity u(x, y, z) on each y plane. We show the physically significant modes with iso-surfaces. As can be seen, the streamwise-elongated streaks and their nonlinear harmonics are extracted by W 1,2 and W 1,3 . W 2,2 and W 3,2 capture streak instabilities at different scales. W 2,3 and W 3,3 stand for the higher-end spectrum in the spanwise direction, displaying the canonical bypass nature of the transition process, i.e., small scales in the freestream only reappear further downstream in the boundary layer leading to a 'clean area'. Comparing the results presented in figure 6 and 8, both inputs of the instantaneous 3D data and its visualisation leads to adaptively and physicallyrelevant decompositions. The extracted streaks and their secondary instabilities from both results match well with each other and correctly reflect the flow physics. Note that when viewed and processed as an image, flow visualisation stands for a projection of the 3D instantaneous data, with the objective to compress while retaining critical features of the data. Differences, as a result, stem from the discrepancy in the Fourier supports. For example, secondary instabilities are extracted by W 2,1 , W 2,2 and W 2,3 from the visualisation, while W 2,2 & W 3,2 from 3D instantaneous data stand for streak instabilities. Another advantage of using flow data is that the extracted modes will have a solid physical meaning. For example, in this case, the modes are the decomposed streamwise velocities. 3D data further leads to more details of the flow, as evidenced by the bypass characteristics in W 2,3 and W 3,3 .

Comparison with the other methods
At this point, it becomes merited to ask, what is the added value of EWT modes to the state-of-art data-based flow decompositions (e.g. the POD and DMD families)? This section takes the boundary layer transitional flow as an example to show the strengths and weaknesses of the proposed methods. The flow case corresponds to the example shown in § 3.2. Time series as inputs for POD and DMD start from the initialisation of FST, encompassing the development of streamwise-elongated streaks, their secondary instabilities, and end at the snapshot shown in figure 8 where hairpin vortices appear. Figure 9 shows the decomposed flows using POD, DMD and EWT, respectively. These modes are illustrated with isosurfaces of the streamwise velocity. From Figure 9(d) it is evident that the first four POD modes cover most of the energy as displayed in panel (a). As can be seen, the most energetic structure is streamwise-elongated streaks. Mode 2 & 3 displayed a mixed structure of streaks and their instabilities. Mode 4, on the other hand, extracts smaller scales The EWT modes are obtained solely from a single snapshot corresponding to figure 5 and 8. Without resorting to time series, the flows are decomposed according to its spatial scales. As a result, the streaks, their secondary instabilities and the bypass nature of smaller scales are distinctly identified.
On the aspect of outputs, all three methods can lead to an arbitrary number of modes. A notable question arises as what is the effectiveness of the modes to characterise the flow. For example, how many modes would be sufficient to represent the flow and what is the criteria for selecting dominant modes? This has been recognised as a weakness of modal based flow analysis (Taira et al. 2017). POD modes are arranged according to their energy, which is a meaningful indication of importance. However, highly nonlinear problems may require a large number of POD modes to cover the majority of the energy (Murata et al. 2020). On the other hand, a single correct way to rank DMD modes is absent though the growth rate and amplitude provide means of measurement. With regard to EWT, there is not a concomitant indicator (e.g. energy, frequency, growth rate) associated with the modes. Instead of producing a large number of modes and choosing a few important ones (e.g. in POD and DMD), EWT aims at generating a limited number of modes, differing only in its spatial scales. We have shown that for flows with compact Fourier supports, M x = M z = 3 gives satisfactory results. However, for flows with a broader spatial spectrum, one may need to increase the number of modes to reduce the mixing of scales in extracted modes.
EWT shares a similar theoretical basis, i.e. the multi-resolution analysis, with CVE (Farge et al. 2001). In CVE, wavelet coefficients are first obtained for the vorticity field using orthogonal wavelet transform. These coefficients are further divided into two groups (coherent and incoherent) according to their amplitudes. The coherent and incoherent flow components are therefore reconstructed through the inverse wavelet transform. By construction, the objective of CVE is to extract coherent structures based on the wavelet denoising method rather than decomposing a flow into several modes. Results of CVE for the boundary layer transitional flow are provided in figure 10. For this case, the coherent structure is captured by only 0.8% of the total wavelet coefficients, reflecting its sparsity distribution in the wavelet space.
As such, it becomes clear that EWT provides a new decomposition strategy, as shown

Input
Decomposition strategy POD flow data as time series by energy DMD flow data as time series by frequency and growth rate CVE instantaneous flow data by amplitude of wavelet coefficients EWT instantaneous flow data or its visualisation by averaged spatial Fourier supports Table 1. A comparison of POD, DMD, CVE and EWT.
in table 1. Furthermore, EWT does not depend on a separation-of-variable strategy, which leads to its advantage dealing with non-stationary flows or traveling wave problems.

Concluding remarks
In this paper, we proposed an image-based flow decomposition using the 2D tensor empirical wavelet transform (EWT) (Gilles 2013;Gilles et al. 2014), which was initially devised for image processing. According to the Fourier supports of the input data, a set of adaptive filter banks (an orthogonal frame) is built to perform the decomposition. A first example considers the interactions between a 2D wake and a jet plume, where only experimental flow visualisations are available. The EWT modes correctly isolate the jet and wake components and their instabilities. In the spirit of 2D tensor EWT, the visualisation's principle direction shall coincide with one of the decomposition directions. We show that this direction can be accurately determined using the shadow mode of EWT (W 1,1 ). The second example considers an early-stage boundary layer transition subject to FST, where direct numerical simulation provided full data-set. For both inputs of 3D instantaneous flow data and its visualisation, the EWT modes distinctly extract streamwise-elongated streaks, their secondary instabilities and smaller scales. A comparison with biGlobal stability analysis justifies the EWT modes that characterise the secondary instabilities. The bypass nature of smaller scales is also captured by EWT modes based on the 3D instantaneous flow data.
Compared to the prevailing data-based methods on flow decomposition (POD, DMD, CVE, to name a few), EWT offers a new strategy to adaptively decompose a flow from its averaged Fourier supports. An instantaneous flow or its visualisation is thus readily decomposed without resorting to its time series. The method also functions well on nonstationary flows or traveling wave problems by extracting fluid physics that are localised from the input. Still, it would be less effective to flows with broader Fourier supports, e.g. fully developed turbulent flows. The number of EWT modes need to be in line with the spatial spectrum of the flow. Future development of the method can be extended to account for temporal spectrum and the prediction of flow evolution.

Declaration of interests
The authors report no conflict of interests.