Multi-Scale Displacement Estimation and Registration for 2-D and - PowerPoint PPT Presentation

Multi-Scale Displacement Estimation and Registration for 2-D and 3-D Datasets Nick Kingsbury Signal Processing Group, Dept. of Engineering University of Cambridge, Cambridge CB2 1PZ, UK. ngk@eng.cam.ac.uk www.eng.cam.ac.uk/~ngk July 2005 UNIVERSITY OF CAMBRIDGE

Displacement Estimation and Registration – 1 Nick Kingsbury The Problem: Efficient Displacement Estimation / Registration of Noisy Data Applications: • Registration of medical datasets taken some time apart and correction for patient movement • Conversion from low-quality video to high-quality still images – e.g. correction of fluctuations in atmospheric refraction (heat shimmer) • Motion estimation for non-rigid objects and fluids • Registration of multi-look images affected by speckle, usually due to illumination from coherent sources such as lasers or synthetic aperture radar (SAR). Displacement estimation usually involves measuring gradients, derivatives or differences . High noise levels mean that registration algorithms must be robust to noise if the noise is uncorrelated between images.

Displacement Estimation and Registration – 2 Nick Kingsbury Key Features of Robust Registration Algorithms • Edge-based methods are more robust than point-based ones. • Bandlimited multiscale (wavelet) methods allow spatially adaptive denoising. • Phase-based bandpass methods can give rapid convergence and immunity to illumination changes between images ( but we have to be careful about 2 π ambiguities) . • If the displacement field is smooth, a wider-area parametric (affine) model of the field is likely to be more robust than a highly-local translation-only model. Note: Biological vision systems have evolved to use multiscale directional bandpass filters as their front-end process (e.g. the V1 cortical filters in humans / mammals).

Displacement Estimation and Registration – 3 Nick Kingsbury Selected Methods • Dual-tree Complex Wavelet Transform (DT CWT): ◦ efficiently synthesises a multiscale directional shift-invariant filterbank, with perfect reconstruction; ◦ provides complex coefficients whose phase shift depends approximately linearly on displacement; ◦ allows each subband of coefficients to be interpolated (shifted) independently of other subbands (because of shift invariance of the subband H ( z ) ). • Parametric model of displacement field, whose solution is based on local edge-based motion constraints (Hemmendorff, Andersson, Kronander and Knutsson, IEEE Trans Medical Imaging, Dec 2002): ◦ derives straight-line constraints from directional subbands of the DT CWT; ◦ solves for spatially-varying affine model parameters which minimise constraint error energy over multiple directions and scales.

Displacement Estimation and Registration – 4 Nick Kingsbury Q-shift Dual Tree Complex Wavelet Transform in 1-D Level 4 ★✥ x 0000 a x 000 a ✲ H 00 a Level 3 ✲ ✲ ↓ 2 Tree a ★✥ ( 3 q ) ✧✦ x 00 a ✲ H 00 a Level 2 ✲ even ↓ 2 ★✥ ★✥ ( 3 q ) ✧✦ x 0001 a x 0 a ✲ H 00 a ✲ H 01 a Level 1 ✲ ✲ ✲ even ↓ 2 ↓ 2 ★✥ ★✥ ( 3 q ) ✧✦ ( q ) ✧✦ ✲ H 0 a ✲ H 01 a ✲ x 001 a ✲ ✲ ↓ 2 even ↓ 2 ★✥ ✘ ✘ (1) ✧✦ ( q ) ✧✦ ✘ ✘ ✘ ✘ ✲ H 01 a ✘ ✘ ✲ x 01 a ★✥ ✲ ✘ ↓ 2 ✘ odd ★✥ ✘ x 0000 b ✘ x 000 b ( q ) ✧✦ ✘ ✘ ✲ H 00 b ✘ ✘ ✲ H 1 a 1-D ✘ ✲ ✲ ✘ ↓ 2 ✲ x 1 a ★✥ ✲ ✘ ↓ 2 ✘ ( q ) ✧✦ ✘ ✘ x 00 b (0) ✧✦ Input ✘ ✘ ✲ H 00 b ✘ ✘ ✘ ✲ ✘ ↓ 2 even ★✥ ★✥ Signal ✘ ✘ ( q ) ✧✦ ✘ x 0001 b ✘ x 0 b ✘ ✘ ✲ H 00 b ✲ H 01 b ✘ ✘ x ✘ ✲ even ✲ ✲ ✘ ↓ 2 ↓ 2 ★✥ ★✥ ✘ ✘ ( q ) ✧✦ ( 3 q ) ✧✦ ✘ ✘ ✲ H 0 b ✲ H 01 b ✲ x 001 b ✲ even ✲ ↓ 2 ↓ 2 ★✥ ✧✦ ✧✦ (0) ( 3 q ) ✲ H 01 b ✲ x 01 b ✲ ↓ 2 odd ★✥ ( 3 q ) ✧✦ Tree b ✲ H 1 b ✲ x 1 b ✲ ↓ 2 (1) ✧✦ Figure 1: Dual tree of real filters for the Q-shift CWT, giving real and imaginary parts of complex coefficients from tree a and tree b respectively. Figures in brackets indicate the approximate delay for each filter, where q = 1 4 sample period.

Displacement Estimation and Registration – 5 Nick Kingsbury Q-shift DT CWT Basis Functions – Levels 1 to 3 1 Level 1: scaling fn. Tree a 0 Tree b wavelet |a + jb| −1 Level 2: scaling fn. −2 wavelet −3 Level 3: scaling fn. −4 wavelet −5 −20 −15 −10 −5 0 5 10 15 20 sample number Figure 2: Basis functions for adjacent sampling points are shown dotted.

Displacement Estimation and Registration – 6 Nick Kingsbury 2-D Basis Functions at level 4 DT CWT real part 0.1 Level 4 Sc. funcs. 0.05 0 15 45 75 −75 −45 −15(deg) −60 −40 −20 0 20 40 60 0.1 Level 4 Wavelets Magnitude 0.05 0 DT CWT imaginary part −0.05 Real Imaginary −60 −40 −20 0 20 40 60 Real DWT Time (samples) e jω 1 x e jω 2 y = e j ( ω 1 x + ω 2 y ) 90 45(?) 0 Figure 3: Basis functions of 2-D Q-shift complex wavelets (top), and of 2-D real wavelet filters (bottom), all illustrated at level 4 of the transforms. The complex wavelets provide 6 directionally selective filters, while real wavelets provide 3 filters, only two of which have a dominant direction. The 1-D bases, from which the 2-D complex bases are derived, are shown to the right.

Displacement Estimation and Registration – 7 Nick Kingsbury Test Image and Colour Palette for Complex Coefficients Colour palette for complex coefs. −1 −0.8 −0.6 −0.4 Imaginary part −0.2 0 0.2 0.4 0.6 0.8 1 −1 −0.5 0 0.5 1 Real part

Displacement Estimation and Registration – 8 Nick Kingsbury 2-D DT-CWT Decomposition into Subbands Figure 4: Four-level DT-CWT decomposition of Lenna into 6 subbands per level (only the central 128 × 128 portion of the image is shown for clarity). A colour-disc palette (see previous slide) is used to display the complex wavelet coefficients.

Displacement Estimation and Registration – 9 Nick Kingsbury Image A Image B Registration ❄ ❄ Algorithm: DT CWT DT CWT ❄ ❄ Select CWT scales (levels) according to iteration ❄ ✲ Shift within subbands t Generate displacement ❄ ❄ field v ( x i ) Form constraints c i ✻ parameter ❄ t field a X Calculate Q X at each locality X Store a X ❄ (initialise to zero) Smooth elements of Q X across image ✻ ✤✜ ❄ ✛ Solve for a X , min at each X + ✣✢ ✻ increment of ❄ parameter field a X Inverse DT CWT ❄ Image A registered to Image B

Displacement Estimation and Registration – 10 Nick Kingsbury Parametric Model: Linear Constraint Equations � � v ( x i ) Let the displacement vector at the i th location x i be v ( x i ) ; and let ˜ v i = . 1 Note that, as well as x i , the locator i also specifies a subband direction d i ( 1 . . . 6 ) and a scale (level) s i . A straight-line constraint on v ( x i ) can be written c T v i = 0 c 1 ,i v 1 ,i + c 2 ,i v 2 ,i + c 3 ,i = 0 or i ˜ For a phase-based system in which wavelet coefficients at { x i , d i , s i } in images A and B have phases θ A,i and θ B,i , approximate linearity of phase θ vs. displacement v ( x i ) means that � � ∇ x θ i c T v i ≈ 0 if c i = C i i ˜ θ B,i − θ A,i In practise we compute this by averaging finite differences at the centre x i of a 2 × 2 × 2 block of coefficients from a given subband { d i , s i } of images A and B . Note: C i is a constant which does not affect the line defined by the constraint, but it is important as a weight for combining constraint errors (see later).

Displacement Estimation and Registration – 11 Nick Kingsbury Parameters of the Model We can define a 6-term affine parametric model a for v such that � a 1 � � a 3 � � x 1 � a 5 v ( x ) = + a 2 a 4 a 6 x 2 or in a more useful form   a 1 � � 1 0 0 0 x 1 x 2 . .  = K ( x ) . a . v ( x ) = .  0 1 0 x 1 0 x 2 a 6 Affine models can synthesise translation, rotation, constant zoom, and shear. A quadratic model , which allows for linearly changing zoom (approx perspective), requires up to 6 additional parameters and columns in K of the form x 2 x 2 � � 0 0 0 . . . x 1 x 2 1 2 x 2 x 2 0 0 0 . . . x 1 x 2 1 2

Displacement Estimation and Registration – 12 Nick Kingsbury Solving for the Model Parameters Using techniques (due to Hemmendorff et al) similar to homogeneous coordinates: � � � � K ( x i ) 0 a Let ˜ v i = ˜ K i = and ˜ a = so that ˜ a . K i ˜ 1 1 0 Ideally for a given scale-space locality X , we wish to find the parametric vector ˜ a such that c T v i = ˜ v i = 0 when for all i such that { x i , d i , s i } ∈ X . i ˜ ˜ K i ˜ a In practise this is an overdetermined set of equations, so we find the LMS solution , i.e. the value of a which minimises the squared error v i || 2 = a || 2 = a T ˜ a T ˜ � � � || c T || c T i ˜ K T i c i c T i ˜ E X = a = ˜ i ˜ K i ˜ ˜ K i ˜ Q X ˜ a i ∈X i ∈X i ∈X � where ˜ K T ˜ i c i c T i ˜ Q X = K i . i ∈X