Synchrosqueezed Curvelet Transforms for 2D mode decomposition - PowerPoint PPT Presentation

Synchrosqueezed Curvelet Transforms for 2D mode decomposition Haizhao Yang Department of Mathematics, Stanford University Collaborators: Lexing Ying † and Jianfeng Lu ♯ † Department of Mathematics and ICME, Stanford University ♯ Department of Mathematics, Duke University SIAM Conference on Imaging Science, May 2014

Geophysics ◮ A superposition of several wave-like components. ◮ Wave field separation and ground roll removal problems. 1 0.8 0.6 0.4 0.2 0 −0.2 −0.4 −0.6 0.4 0.5 0.3 0.4 0.2 0.3 0.1 0.2 0 0.1 −0.1 0 −0.2 −0.1 −0.3 −0.2 −0.4 −0.3

Materials science Atomic crystal analysis ◮ Observation: an assemblage of wave-like components; ◮ Goal: Crystal segmentations, crystal rotations, crystal defects, crystal deformations. Figure : Left: An atomic crystal image. Middle: Estimated crystal rotation. Right: Identified crystal defects. Takes about 10s for a 1024*1024 image, while other methods using GPU computing take at least 40s.

1D mode decomposition Known: A superposition of wave-like components � K � K α k ( t ) e 2 π iN k φ k ( t ) . f ( t ) = f k ( t ) = k =1 k =1 Unknown: Number K , components f k ( t ), smooth instantaneous amplitudes α k ( t ), smooth instantaneous frequencies N k φ ′ k ( t ). Existing methods: ◮ Empirical mode decomposition methods (Huang et al. 98, 09); ◮ Synchrosqueezed wavelet transform (Daubechies et al. 09, 11); Synchrosqueezed wave packet transform (Y. 13); ◮ Data-driven time-frequency analysis (Hou et al. 11, 12, 13); ◮ Regularized nonstationary autoregression (Fomel 13);

1D synchrosqueezed wavelet transform (SSWT) Continuous wavelet transform of a signal s ( t ): � s ( t ) a − 1 / 2 φ ( t − b W s ( a , b ) = � s ( t ) , φ ab ( t ) � = ) d t . a EX: s ( t ) = A cos( ω t ). � 1 φ ( a ξ ) e ib ξ d ξ = A s ( ξ ) a 1 / 2 � 4 π a 1 / 2 � φ ( a ω ) e ib ω . W s ( a , b ) = � 2 π Synchrosqueezing for better readability Figure : Numerical examples by Daubechies et al, signal f ( t ) = sin (8 t ).

Definitions of SSWT EX: s ( t ) = A cos( ω t ). 4 π a 1 / 2 � φ ( a ω ) e ib ω ⇒ ∂ b W s ( a , b ) A W s ( a , b ) = iW s ( a , b ) = ω Definition: Instantaneous frequency estimate ω s ( a , b ) = ∂ b W s ( a , b ) iW s ( a , b ) . Definition: Synchrosqueezed wavelet transform � W s ( a , b ) a − 3 / 2 δ ( ω s ( a , b ) − ω ) d a . T s ( ω, b ) = { a : W s ( a , b ) � =0 }

Theory of SSWT Theorem: (Daubechies, Lu, Wu 11 ACHA) If � K � K α k ( t ) e 2 π iN k φ k ( t ) f ( t ) = f k ( t ) = k =1 k =1 and f k ( t ) are well-separated, then ◮ T f ( a , b ) has well-separated supports Z k concentrating ( N k φ ′ k ( b ) , b ); ◮ f k ( t ) can be accurately recovered by applying an inverse transform on I Z k ( a , b ) T f ( a , b ). where I Z k ( a , b ) is an indication function.

2D mode decomposition Known: A superposition of wave-like components K K � � α k ( x ) e 2 π iN k φ k ( x ) . f ( x ) = f k ( x ) = k =1 k =1 Unknown: Number K , components f k ( x ), local amplitudes α k ( x ), local wave vectors N k ∇ φ k ( x ). Existing methods: 1. 2D EMD methods (Huang et al.); 2. 2D Empirical wavelet, curvelet transforms (Gilles, Tran, Osher 13); 3. 2D SS wave packet and curvelet transforms (Y. and Ying 13, 14).

2D wavelet transform or wave packet transform? Consider two plane waves e ip · x and e iq · x again. 10 10 10 8 8 8 6 6 6 4 4 4 2 2 2 ξ 2 0 ξ 2 0 ξ 2 0 −2 −2 −2 −4 −4 −4 −6 −6 −6 −8 −8 −8 −10 −10 −10 −10 −5 0 5 10 −10 −5 0 5 10 −10 −5 0 5 10 ξ 1 ξ 1 ξ 1 ◮ Left: Supports of continuous wavelets and plane waves with different wave numbers in the Fourier domain. Radial separation. ◮ Middle: Supports of wavelets and plane waves with the same wave number. No angular separation. ◮ Right: Supports of wave packets and plane waves with the same wave number. Angular separation.

Definition: Wave Packets Given the mother wave packet w ( x ) and s ∈ (1 / 2 , 1), the family of wave packets { w pb ( x ) , p , b ∈ R 2 } is defined as w pb ( x ) = | p | s w ( | p | s ( x − b )) e 2 π i ( x − b ) · p , or equivalently in Fourier domain w pb ( ξ ) = | p | − s e − 2 π ib · ξ ˆ w ( | p | − s ( ξ − p )) . � Definition: The 2D wave packet transform of a function f ( x ) is a function of p , b ∈ R 2 � W f ( p , b ) = w pb ( x ) f ( x ) d x .

Definition: The local wavevector estimation of a function f ( x ) at ( p , b ) is v f ( p , b ) = ∇ b W f ( p , b ) 2 π iW f ( p , b ) for p , b ∈ R 2 such that W f ( p , b ) � = 0. Definition: Given f ( x ), for v , b ∈ R 2 , W f ( p , b ), and v f ( p , b ), the synchrosqueezed energy distribution T f ( v , b ) is � | W f ( p , b ) | 2 δ ( v f ( p , b ) − v ) d p . T f ( v , b ) = { p : W f ( p , b ) � =0 } Remark: The support of T f ( v , b ) is concentrating around local wavevectors.

Sketch of 2D mode decomposition Two components: f ( x ) = α 1 ( x ) e 2 π iN φ 1 ( x ) + α 2 ( x ) e 2 π iN φ 2 ( x ) . 200 200 200 200 100 100 100 100 p 2 0 v 2 0 v 2 0 v 2 0 −100 −100 −100 −100 −200 −200 −200 −200 200 200 200 200 100 1 100 1 100 1 100 1 0 0 0 0 −100 0.5 −100 0.5 −100 0.5 −100 0.5 −200 −200 −200 −200 p 1 0 v 1 0 v 1 0 v 1 0 b 2 b 2 b 2 b 2 Fast algorithms for forward and inverse transforms. O ( L 2 log L ) for a L × L image.

Theory of 2D SS wave packet transform (SSWPT) Theorem [Y.,Ying 13, SIAM Imaging Science] For a well-separated superposition f ( x ) = � K k =1 α k ( x ) e 2 π iN φ k ( x ) and ǫ > 0, we define R f ,ǫ = { ( p , b ) : | W f ( p , b ) | ≥ | p | − s √ ǫ } and Z f , k = { ( p , b ) : | p − N ∇ φ k ( b ) | ≤ | p | s } for 1 ≤ k ≤ K . There exists a constant ǫ 0 ( K ) > 0 such that for any ǫ ∈ (0 , ǫ 0 ) there exists a constant N 0 ( K , ǫ ) > 0 such that for any N > N 0 ( K , ǫ ) the following statements hold. (i) { Z f , k : 1 ≤ k ≤ K } are disjoint and R f ,ǫ ⊂ � 1 ≤ k ≤ K Z f , k ; (ii) For any ( p , b ) ∈ R f ,ǫ ∩ Z f , k , � √ ǫ. | v f ( p , b ) − N ∇ φ k ( b ) | | N ∇ φ k ( b ) |

Banded wave-like components 1 1 12 0.9 0.8 0.9 10 0.8 0.6 0.8 0.7 0.4 0.7 8 0.2 0.6 0.6 x 2 0.5 0 x 2 0.5 6 0.4 −0.2 0.4 4 0.3 −0.4 0.3 0.2 −0.6 0.2 2 0.1 −0.8 0.1 0 −1 0 0 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x 1 x 1 1 1 0.07 0.6 0.9 0.9 0.06 0.8 0.8 0.5 0.7 0.7 0.05 0.4 0.6 0.6 0.04 x 2 0.5 x 2 0.5 0.3 0.03 0.4 0.4 0.2 0.3 0.3 0.02 0.2 0.2 0.1 0.01 0.1 0.1 0 0 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x 1 x 1 Figure : Top-left: A banded deformed plane wave. Top-right: Number of wave packet coefficients | W f ( a , b ) | > δ . Bottom-left: Relative error of local wave-vector estimates using SSWPT. Bottom-right: Relative error of local wave-vector estimates using SSCT.

General curvelet transform: Some notations 1. The scaling matrix � a t � 0 A a = . a s 0 2. The rotation angle θ and rotation matrix � cos θ � − sin θ R θ = . sin θ cos θ 3. The unit vector e θ = (cos θ, sin θ ) T of rotation angle θ . 4. θ α represents the argument of given vector α . Definition: 2D general curvelets For 1 2 < s < t < 1, define t + s 2 e 2 π ia ( x − b ) · e θ w ( A a R − 1 w a θ b ( x ) = a θ ( x − b )) . A family of curvelets { w a θ b ( x ) , a ∈ [1 , ∞ ) , θ ∈ [0 , 2 π ) , b ∈ R 2 } .

Definition: 2D general curvelet transform � W f ( a , θ, b ) = R 2 w a θ b ( x ) f ( x ) d x for a ∈ [1 , ∞ ), θ ∈ [0 , 2 π ), b ∈ R 2 . Definition: local wave vector estimation Local wave-vector estimation at ( a , θ, b ) is v f ( a , θ, b ) = ∇ b W f ( a , θ, b ) 2 π iW f ( a , θ, b ) for a ∈ [1 , ∞ ), θ ∈ [0 , 2 π ), b ∈ R 2 such that W f ( a , θ, b ) � = 0. Definition: synchrosqueezed energy distribution Synchrosqueezed energy distribution T f ( v , b ) is � | W f ( a , θ, b ) | 2 δ ( ℜ v f ( a , θ, b ) − v ) a d a d θ T f ( v , b ) = { ( a ,θ ): W f ( a ,θ, b ) � =0 } for v ∈ R 2 , b ∈ R 2 .

Theory of the SS curvelet transform (SSCT) Theorem [Y., Ying SIAM Mathematical Analysis 14] Suppose 1 2 < s < η < t are fixed, and a well-separated superposition � K e − ( φ k ( x ) − c k ) 2 /σ 2 k α k ( x ) e 2 π iN φ k ( x ) , f ( x ) = k =1 where σ k ≥ N − η . For any ǫ > 0, define � � 2 √ ǫ ( a , θ, b ) : | W f ( a , θ, b ) | ≥ a − s + t R f ,ǫ = and � � a R − 1 ( a , θ, b ) : | A − 1 Z f , k = θ ( a · e θ − N ∇ φ k ( b )) | ≤ 1 for 1 ≤ k ≤ K . For any ǫ , there exists N 0 ( ǫ ) > 0 such that for any N > N 0 ( ǫ ) the following statements hold. (i) { Z f , k : 1 ≤ k ≤ K } are disjoint and R f ,ǫ ⊂ � 1 ≤ k ≤ K Z f , k ; (ii) For any ( a , θ, b ) ∈ R f ,ǫ ∩ Z f , k , � √ ǫ. | v f ( a , θ, b ) − N ∇ φ k ( b ) | | N ∇ φ k ( b ) |

Applications in Geophysics Noisy Example 1 for SS curvelet transform 1 0.5 0 −0.5 −1 0.6 0.4 0.5 0.3 0.4 0.2 0.3 0.1 0.2 0 0.1 −0.1 0 −0.2 −0.1 −0.3 −0.2 −0.4 −0.3 Figure : Top: A noisy superposition of two components. Left: First recovered component. Right: Second recovered component.