COMPRESSIVE OPTICAL DEFLECTOMETRIC TOMOGRAPHY Adriana Gonz alez - PowerPoint PPT Presentation

COMPRESSIVE OPTICAL DEFLECTOMETRIC TOMOGRAPHY Adriana Gonz´ alez ICTEAM/UCL March 26th, 2014 1

ISPGroup - ICTEAM - UCL Universit´ e catholique de Louvain, Louvain-la-Neuve, Belgium. ISP Group 4 Professors 17 researchers http://sites.uclouvain.be/ispgroup 12 PhD students Compressed Sensing Group Prof. Laurent Prof. Christophe De Dr. Prasad K´ evin Degraux Jacques Vleeschouwer Sudhakar 2

Outline Optical Deflectometric Tomography 1 Compressiveness in RIM Reconstruction 2 Compressiveness in Acquisition 3 3

Outline Optical Deflectometric Tomography 1 Compressiveness in RIM Reconstruction 2 Compressiveness in Acquisition 3 4

Optical Deflectometric Tomography Interest Optical characterization of (transparent) objects ODT Tomographic Imaging Modality Measures light deviation caused by the difference in the object refractive index e 2 Deviated Light Rays α ( τ , θ ) p θ t θ θ e 1 τ n ( r ) Incident Light Rays 5

Schlieren Deflectometer SLM Rotating CCD Lens 1 ( f ) Lens 2 Lens 3 (modulation by ) ϕ i object Intensity p Pinhole change α Uniform ∆ x = f tan α τ Light e 2 Optical axis Source e 3 − τ I 0 ∆ x Telecentric − θ y p system e 1 s p y p = � s p , ϕ i � 6

Schlieren Deflectometer s p y p = � s p , ϕ i � 1 Compressiveness in ϕ sinusoidal pattern ⇒ 4 shifted patterns ϕ 1 , ϕ 2 , ϕ 3 , ϕ 4 ⇒ 4 measurements to recover α RIM reconstruction Assuming deflections at one point Objects RIM variation only on e 1 − e 2 ⇒ α , 2-D slices 7

Schlieren Deflectometer s p y p = � s p , ϕ i � 2 Compressiveness in Deflections produced by several points acquisition Objects RIM variation also on e 3 ⇒ α and β , 3-D volume M binary modulation patterns ϕ i to eliminate nonlinearities 8

Outline Optical Deflectometric Tomography 1 Compressiveness in RIM Reconstruction 2 Compressiveness in Acquisition 3 9

Framework Joint work with Prof. Laurent Jacques and Prof. Christophe De Vleeschouwer from UCL and Dr. Philippe Antoine from Lambda-X Problem To reconstruct the refractive index map of transparent materials from light deflection measurements ( α ) under few orientations ( θ ) Assumption Objects are constant along the e 3 direction Deflections at only one point [1] A. Gonz´ alez et al. iTWIST12 [2] P. Antoine et al. OPTIMESS 2012 [3] A. Gonz´ alez et al. IPI Journal (2014) 10

Continuous facts Mathematical Model e 2 Deviated Light Rays Eikonal equation � � α ( τ, θ ) d n d R curved : r ( s ) → d s r ( s ) = ∇ n p θ d s t θ Approximation θ small α → R straight : r · p θ = τ e 1 τ error < 10% n ( r ) ∆( τ, θ ) = sin( α ) Incident Light Rays � � � 1 δ ( τ − r · p θ ) d 2 r ∆( τ, θ ) = ∇ n ( r ) · p θ R 2 n r Deflectometric Central Slice Theorem � � � 2 π i ω R ∆( τ, θ ) e − 2 π i τ ω d τ = y ( ω, θ ) := n r � ω p θ n � � � : 2-D Fourier transform of � ω p θ n in Polar grid n 11

Discrete Forward Model 2 π i ( δ r ) 2 diag( ω (1) , · · · , ω ( M ) ) � y = n n r ⇓ y = DF n + η n ∈ R N ; Cartesian grid of N = N 2 0 pixels; sampling: δ r y ∈ R M ; Polar grid of M = N τ N θ pixels; sampling: δτ , δθ D : 2 π i ( δ r ) 2 diag( ω (1) , · · · , ω ( M ) ) ∈ C M × M n r F ∈ C M × N : Non-equispaced Fourier Transform (NFFT) [4] η ∈ C M : numerical computations, model discretization, model discrepancy, observation noise [4] J. Keiner et al. (2009) 12

ODT vs. AT y = DF n + η Main difference: Operator D Without noise η → classical tomographic model y = D − 1 y = F n ˜ For η � = 0 → Not a classical tomographic model - η : AWGN → D − 1 η not homoscedastic 13

ODT vs. AT Observation: 1-D FT of sinograms along the τ direction 0.12 −3 0.1 0.1 −2 0.08 0.08 0.06 −1 Absorption 0.04 0.06 y( ω , θ ) ω 0 0.02 Tomography 0.04 0 1 0.02 −0.02 2 0 −0.04 3 −0.06 −0.02 0 50 100 150 −4 −3 −2 −1 0 1 2 3 4 θ ω −3 −3 4 x 10 x 10 −3 8 3 6 −2 Optical 2 4 −1 Deflection 2 y( ω , θ ) 1 ω 0 Tomography 0 0 1 −2 −1 2 −4 −2 3 0 50 100 150 θ −4 −3 −2 −1 0 1 2 3 4 ω 14

Naive Reconstruction Methods y = Φ n + η = DF n + η Filtered Back Projection Analytical method Solution ˜ n FBP : - Filtering the tomographic projections AT: ramp filter; ODT: Hilbert filter - Backprojecting in the spatial domain by angular integration Minimum Energy Reconstruction n ME = Φ † y = Φ ∗ ( ΦΦ ∗ ) − 1 y ˜ ˜ ≡ n ME = arg min � u � 2 s . t . y = Φu u ∈ R N Noise Solution: Problems: Compressiveness ⇒ M ( N θ ) < N Regularization ⇒ ill-posed problem 15

Sparsity prior Heterogeneous transparent materials with slowly varying refractive index separated by sharp interfaces TV and BV promote the perfect “cartoon shape” model “Sparse” gradient ⇓ Small Total Variation norm � n � TV := � ∇ n � 2 , 1 16

Other priors Positive RIM ⇒ n � 0 The object is completely contained in the image. Pixels in the border are set to zero in order to guarantee uniqueness of the solution. ⇒ n | δ Ω = 0 SOLUTION UNIQUENESS 17

TV- ℓ 2 reconstruction and Noise y = Φ n + η = DF n + η TV- ℓ 2 Reconstruction n TV − ℓ 2 = arg min ˜ � u � TV s . t . � y − Φu � 2 ≤ ε, u � 0 , u ∂ Ω = 0 u ∈ R N Noise Observation noise → σ 2 obs Modeling error → ray tracing with Snell law ≈ 10% Interpolation noise → NFFT error (very small) 18

TV- ℓ 2 reconstruction y = Φ n + η = DF n + η TV- ℓ 2 Reconstruction ˜ n TV − ℓ 2 = arg min � u � TV s . t . � y − Φu � 2 ≤ ε, u � 0 , u ∂ Ω = 0 u ∈ R N ˜ n TV − ℓ 2 = arg min � u � TV + ı C ( Φu ) + ı P 0 ( u ) u ∈ R N Indicator function: ı X ( x ) = 0 if x ∈ X ; + ∞ otherwise ı C and ı P 0 are the indicator functions into the following convex sets: - C = { v ∈ C M : � y − v � ≤ ε } - P 0 = { u ∈ R N : u i ≥ 0 if i ∈ int Ω; u i = 0 if i ∈ ∂ Ω } Sum of 3 proper, lower semicontinuous, convex functions Reconstruction using CP algorithm [5] expanded in a product space [5] A. Chambolle and T. Pock. Journal of Mathematical Imaging and Vision. (2011) 19

Reconstruction Algorithm Chambolle-Pock (CP)  v ( k +1) = prox σ F ⋆ ( v ( k ) + σ K ¯ x ( k ) ) min x ∈X F ( Kx ) + G ( x )   x ( k +1) = prox τ G ( x ( k ) − τ K ∗ v ( k +1) ) | 2 < 1  F , G : proper, lsc, convex; τσ | | | K | |  x ( k +1) = 2 x ( k +1) − x ( k ) ¯ Proximal mapping z ) + 1 z − z � 2 f : proper, lsc, convex ⇒ prox σ f z = arg min ¯ z σ f (¯ 2 � ¯ 2 e.g., prox σℓ 1 z = SoftTh( z , σ ) Conjugate function F ⋆ ( v ) = max ¯ v � v , ¯ v � − F (¯ v ) CP Product-Space Expansion t �  � x ( k ) � min F j ( K j x ) + G ( x ) v ( k +1) v ( k )  = prox σ F ⋆ + σ K j ¯ , j ∈ { 1 , · · · t }   j j x ∈X � t j j =1 x ( k +1) = prox τ t H ( x ( k ) − τ i v ( k +1) j =1 K ∗ ) j  t   x ( k +1) = 2 x ( k +1) − x ( k ) K = diag( K 1 , · · · , K t ) ¯ 20

Results: TV − ℓ 2 vs. ME & FBP Compressiveness and noise robustness 45 100 40 35 80 30 60 25 RSNR [dB] RSNR [dB] 20 40 15 10 20 TV−L2 20dB 5 TV−L2 10dB TV−L2 ME 20dB 0 0 ME 10dB ME FBP 20dB −5 FBP FBP 10dB −20 −10 0 20 40 60 80 100 0 20 40 60 80 100 N θ / 360 [%] N θ / 360 [%] � ∆ � 2 � n � 2 MSNR = 20 log 10 RSNR = 20 log 10 � η � 2 � n − ˜ n � 2 21

Results: TV- ℓ 2 vs. ME 0.012 0.01 No measurement noise (MSNR = ∞ ) 0.008 0.006 N θ / 360 = 25% 0.004 0.002 0 −3 x 10 0.012 12 0.01 10 0.008 8 6 0.006 4 0.004 2 0.002 0 0 ˜ ˜ n TV − ℓ 2 : RSNR = 71dB n ME : RSNR = 13dB 22

Results: TV- ℓ 2 vs. ME 0.012 0.01 No measurement noise (MSNR = ∞ ) 0.008 0.006 N θ / 360 = 5% 0.004 0.002 0 −3 x 10 0.012 10 0.01 8 0.008 6 0.006 4 2 0.004 0 0.002 −2 0 ˜ ˜ n TV − ℓ 2 : RSNR = 67dB n ME : RSNR = 5dB 23

Results: TV- ℓ 2 vs. ME 0.012 0.01 MSNR = 20dB 0.008 0.006 N θ / 360 = 25% 0.004 0.002 0 −3 −3 x 10 x 10 12 12 10 10 8 8 6 6 4 4 2 2 0 0 ˜ ˜ n TV − ℓ 2 : RSNR = 39dB n ME : RSNR = 13dB 24

Results: TV − ℓ 2 Convergence 0 . 9 Non-Adaptive: step-sizes constant along the iterations → σ = τ = | | | K | | | Adaptive: step-sizes σ and τ are updated according to the residuals of the algorithm [6] −1 10 50 Adapt Non−Adapt −2 45 10 40 −3 10 ||x (k+1) − x (k) || / || x (k) || 35 RSNR [dB] −4 10 30 −5 10 25 −6 10 20 −7 10 15 Adapt Non−Adapt −8 10 10 0 1000 2000 3000 4000 5000 6000 7000 0 1000 2000 3000 4000 5000 6000 7000 # iter # iter [6] T. Goldstein et al. Preprint (2013) 25

Experimental Results 0.06 50 0.05 100 Bundle of 10 fibers immersed in an optical fluid 0.04 150 0.03 200 0.02 Slices MSNR ≈ 10dB 250 0.01 300 0 350 N θ = 60 ⇒ N θ / 360 = 17% −0.01 400 −0.02 450 −0.03 500 −0.04 100 200 300 400 500 600 τ −3 −3 x 10 14 x 10 9 TV−L2 3 Expected 8 12 7 2.5 10 6 8 2 y (mm) 5 6 1.5 4 4 3 1 2 2 0.5 0 1 0 0 −2 0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3 3.5 x (mm) x (mm) 26