Image stitching

Image stitching

Digital Visual Effectsg Yung-Yu Chuang

with slides by Richard Szeliski, Steve Seitz, Matthew Brown and Vaclav Hlavac

Image stitching

• Stitching = alignment + blending

geometrical photometric registration registration

Applications of image stitching

• Video stabilization

Vid i i

• Video summarization

• Video compression

• Video matting

• Panorama creation

• Panorama creation

Video summarization

Video compression

Object removal

input video

Object removal

remove foreground

Object removal

estimate background

Object removal

background estimation

Panorama creation

Why panorama?

• Are you getting the whole picture?

C t C FOV 50 35°

– Compact Camera FOV = 50 x 35°

Why panorama?

• Are you getting the whole picture?

C t C FOV 50 35°

– Compact Camera FOV = 50 x 35°

– Human FOV = 200 x 135°

Why panorama?

• Are you getting the whole picture?

C t C FOV 50 35°

– Compact Camera FOV = 50 x 35°

– Human FOV = 200 x 135°

P i M i 360 180°

– Panoramic Mosaic = 360 x 180°

Panorama examples

• Like HDR, it is a topic of computational

photography seeking ways to build a better photography, seeking ways to build a better camera mostly in software.

M t h d

• Most consumer cameras have a panorama mode

• Mars:

http://www.panoramas.dk/fullscreen3/f2_mars97.html

• Earth:

http://www.panoramas.dk/new-year-2006/taipei.html http://www.360cities.net/

http://www.360cities.net/

What can be globally aligned?

• In image stitching, we seek for a matrix to

globally warp one image into another Are any globally warp one image into another. Are any two images of the same scene can be aligned this way?

this way?

– Images captured with the same center of j ti

projection

– A planar scene or far-away scene

A pencil of rays contains all views

real synthetic

camera y

camera

Can generate any synthetic camera viewg y y

as long as it has the same center of projection!

Mosaic as an image reprojection

mosaic projection plane

• The images are reprojected onto a common plane

• The mosaic is formed on this planeThe mosaic is formed on this plane

• Mosaic is a synthetic wide-angle camera

Changing camera center

• Does it still work? synthetic PP PP1

PP2

Planar scene (or a faraway one)

PP3 PP1

PP2

• PP3 is a projection plane of both centers of projection, so we are OK!

• This is how big aerial photographs are made

Motion models

• Parametric models as the assumptions on the relation between two images

relation between two images.

2D Motion models

Motion models

Translation Affine Perspective 3D rotation

2 unknowns 6 unknowns 8 unknowns 3 unknowns

A case study: cylindrical panorama

• What if you want a 360 field of view?

mosaic projection cylinder

Cylindrical panoramas

• Steps

– Reproject each image onto a cylinder – Blend

– Output the resulting mosaic

Cylindrical panorama

1. Take pictures on a tripod (or handheld)

2 W li d i l di

2. Warp to cylindrical coordinate 3. Compute pairwise alignments

4. Fix up the end-to-end alignment 5 Blending

5. Blending

6. Crop the result and import into a viewer

It is required to do radial distortion correction for better stitching results!

correction for better stitching results!

Taking pictures

K id i t i d h d

Kaidan panoramic tripod head

Translation model

Try to align this in PaintShop Pro

Where should the synthetic camera be

real synthetic

camera y

camera

• The projection plan of some camera

• Onto a cylinder

• Onto a cylinder

Cylindrical projection

Adopted from http://www.cambridgeincolour.com/tutorials/image-projections.htm

Cylindrical projection

Cylindrical projection

Adopted from http://www.cambridgeincolour.com/tutorials/image-projections.htm

Cylindrical projection

y

unwrapped cylinder

x

y

θ

x

f

θ

Cylindrical projection

y

unwrapped cylinder

x

y

z f

θ

f

x

Cylindrical projection

y

unwrapped cylinder

x

y

z f

s=f

f

x

Cylindrical reprojection

top-down view^p Focal lengthFocal length – the dirty secret…– the dirty secret

f = 180 (pixels) f = 280 f = 380 Image 384x300

A simple method for estimating f

w

d f

p

f

Or, you can use other software, such as AutoStich, to help.

to help.

Input images

Cylindrical warping

Blending

• Why blending: parallax, lens distortion, scene motion exposure difference

motion, exposure difference

Blending

Blending

Blending

Assembling the panorama

• Stitch pairs together blend then crop

• Stitch pairs together, blend, then crop

Problem: Drift

• Error accumulation

– small errors accumulate over time

Problem: Drift

(x₁,y₁)

(x_n,y_n)

• Solution ^• copy of first

• Solution

– add another copy of first image at the end

– there are a bunch of ways to solve this problem

image

– there are a bunch of ways to solve this problem

• add displacement of (y₁ – y_n)/(n -1) to each image after the first

• compute a global warp: y’ = y + ax

• run a big optimization problem, incorporating this constraint

– best solution, but more complicated – known as “bundle adjustment”

End-to-end alignment and crop

Viewer: panorama

+ +

+ +

+ + +

+

+ +

example: http://www.cs.washington.edu/education/courses/cse590ss/01wi/projects/project1/students/dougz/index.html

Viewer: texture mapped model

example: http://www.panoramas.dk/

Cylindrical panorama

1. Take pictures on a tripod (or handheld)

2 W li d i l di

2. Warp to cylindrical coordinate 3. Compute pairwise alignments

4. Fix up the end-to-end alignment 5 Blending

5. Blending

6. Crop the result and import into a viewer

Determine pairwise alignment?

• Feature-based methods: only use feature points to estimate parameters

to estimate parameters

• We will study the “Recognising panorama”

paper published in ICCV 2003

• Run SIFT (or other feature algorithms) for each

• Run SIFT (or other feature algorithms) for each image, find feature matches.

Determine pairwise alignment

• p’=Mp, where M is a transformation matrix, p and p’ are feature matches

and p are feature matches

• It is possible to use more complicated models

h ffi ti

such as affine or perspective

• For example, assume M is a 2x2 matrix



 







 



 



 





y x m

m

m m

y x

22 21

12 11

' '

• Find M with the least square error

 

n 2











 y m₂₁ m₂₂ y

 





i

p Mp

1

' 2

Determine pairwise alignment











 x' m₁₁ m₁₂ x x₁m₁₁  y₁m₁₂  x₁^'









 

 

 y' m₂₁ m₂₂ y ^'

1 22

1 21

1m y m y

x  

• Overdetermined system



 

 0 0 ^'





















' 1 1 11

1 1

1 1

0 0

0 0

y x m

y x

y x







 















 ^'

12 2 2

2 0 0 x

m m y

x

 



























' 22

21

0

0 _n

n

n m x

m y

x

 













 0 0 x_n y_n y_n^'

Normal equation

Given an overdetermined system

b Ax 

the normal equation is that which minimizes the sum of the square differences between left and sum of the square differences between left and right sides

b A

Ax

A



Why?

Normal equation

 

)

( x  Ax  b

E ^{( )}  







 a

a

b

 



 



 

 

 



 



b

x a

a

: :

:

...

1 1

11

 



 

 

 



 



 



 



: : :

:

1

 



 

 



 

 



 

 

 

 



m

b

x :

: :







 a

... a

b

nxm n equations m variables

nxm, n equations, m variables

Normal equation





 

 











































 

 j

j j j

j

j a x b

x b

a ¹

1 1 1

1 1

: : :







  



 



 







 













 

j

j ij i

m

j

j

ij x b a x b

a

1 1

: b

Ax





















 



 









 

m j

m j

b

1 1

: : :





 



  



 



 

 



 





 

  n

m

j

j nj n

j

j

nj a x b

x b a

1 1

 ^  ^   ^{ } _ ^ _ ^{   }



a

x

b

E

2

)

( x  Ax b   

 

 

 

  

 



i

i j

j

ij

x b

a E

1 1

)

( x Ax b

Normal equation

   

^{   }

^{  b  }

E

2

)

(  A b   

 

 

 

   

 

 





i

i j

j

ij

x b

a E

1 1

)

( x Ax b



  

 

 

0 x

E

2

n

i m

j

ij

x b a

  a

 





 



  



 





 x



 

1 1

i j

j

 





 

 





 





i

i i i

j j

ij

i

a x a b

a

1

1

1 1

1

2

2

) (

2

0   A

Ax  A

b

 E

b A

Ax

A



→ )

( 2

0 A Ax A b

x 

  ^→ A Ax  A b

Normal equation

 ^{Ax  b} 

   

   ^{Ax  b Ax}  ^{  b }



T T

T

    ^{ }

 ^{x Ax}

^A

^{b b}

 ^{ Ax Ax b b }









  ^{ }

b b b

A x

Ax b

Ax A

x

b Ax

b A

x

T T

T T

T

T

  









    ^{A b x A b x b b}

Ax A

x









b A Ax

A

2

2 

 

E

 x

Determine pairwise alignment

• p’=Mp, where M is a transformation matrix, p and p’ are feature matches

and p are feature matches

• For translation model, it is easier.

 ^{   } 

 













 ⁿ

i

i i

i

i x m y y

x m

E

1

' 2 2

' 2 1

0  E

1

0  m

• What if the match is false? Avoid impact of outliers.

RANSAC

• RANSAC = Random Sample Consensus

A l i h f b fi i f d l i h

• An algorithm for robust fitting of models in the presence of many data outliers

• Compare to robust statistics

• Given N data points x_i, assume that mjority of them are generated from a model with

them are generated from a model with parameters , try to recover .

RANSAC algorithm

Run k times:

(1) d l d l

How many times?

How big?

(1) draw n samples randomly

(2) fit parameters  with these n samples

How big?

Smaller is better

(3) for each of other N-n points, calculate

its distance to the fitted model count the its distance to the fitted model, count the number of inlier points, c

O t t  ith th l t Output  with the largest c

How to define?

Depends on the problem.

How to determine k

p: probability of real inliers

P: probability of success after k trials P: probability of success after k trials

k

p

P  1  ( 1  p ) P 1 ( 1 )

n samples are all inliers a failure

f il ft k t i l failure after k trials

) 1

l ( P

n p k

) 1

log(

) 1

log(

p

k P



 

6 0.6 97

for P=0.99

) 1

log( p

6 0.5 293

Example: line fitting

Example: line fitting

n=2 n=2

Model fitting

Measure distances

Count inliers

c=3 c=3

Another trial

c=3 c=3

The best model

c=15 c=15

RANSAC for Homography

RANSAC for Homography

RANSAC for Homography

Applications of panorama in VFX

• Background plates I b d li h i

• Image-based lighting

Troy (image-based lighting)

http://www.cgnetworks.com/story custom.php?story id=2195&page=4 http://www.cgnetworks.com/story_custom.php?story_id 2195&page 4

Spiderman 2 (background plate)

Cylindrical projection

– Map 3D point (X,Y,Z) onto a unit cylinder

unit cylinder

Y

Z X – Convert to cylindrical coordinates

unit cylinder

– Convert to cylindrical image di t

coordinates

unwrapped cylinder

cylindrical image

3D → 2D perspective projection

(X_c,Y_c,Z_c) ( _c, _c, _c) u_c

f

u

Reference

• Richard Szeliski, Image Alignment and Stitching, unpublished draft, 2005.

• R. Szeliski and H.-Y. Shum. Creating full view panoramic image mosaics and texture-mapped models, SIGGRAPH 1997, pp251-258.

M Brown D G Lowe Recognising Panoramas ICCV 2003

• M. Brown, D. G. Lowe, Recognising Panoramas, ICCV 2003.

Direct vs feature-based

• Direct methods use all information and can be very accurate but they depend on the fragile very accurate, but they depend on the fragile

“brightness constancy” assumption

• Iterative approaches require initialization

• Iterative approaches require initialization

• Not robust to illumination change and noise images

images

• In early days, direct method is better.

• Feature based methods are now more robust and potentially faster

and potentially faster

• Even better, it can recognize panorama without initialization

initialization

TODO

• Bundle adjustment LM h d

• LM method

• Direct method vs feature-based method

• Frame-rate image alignment for stabilization

• Frame rate image alignment for stabilization Ri k’ CGA 1995 ? LM th d

• Rick’s CGA 1995 paper? LM method

Project #2 Image stitching

• camera availability T i d?

• Tripod?

• http://www.tawbaware.com/maxlyons/

• http://www.cs.washington.edu/education/cou rses/cse590ss/CurrentQtr/projects.htmQ p j

• http://www.cs.ubc.ca/~mbrown/panorama/pa norama html

norama.html

blending

• Alpha-blending Ph

• Photomontage

• Poisson blending

• Adelson’s pyramid blending

• Hdr?

• Hdr?

3D interpretation



 







 



 

1 0 1 0

0 R t p u K

p

































0 0 0

0 0 0

1 0 1 0

23 22

21

13 12

11

Y X t

r r

r

t r

r r

f

f _x





 







 







 







1 1

0 0

0 1

0 0 0

0 1 0 0

0 0 0

33 32

31

23 22

21

Z Y t

r r

r

t r

r r

f

z y













u

Cylindrical warping

• Given focal length f and image center (x_c,y_c)

(X Y Z)

Y

(X,Y,Z)

(sin,h,cos)

X Y

Z

( )

Cylindrical projection

– Map 3D point (X,Y,Z) onto cylinder

cylinder

Y

Z X – Convert to cylindrical coordinates

unit cylinder

– Convert to cylindrical image di t

coordinates

unwrapped cylinder

cylindrical image

Cylindrical reprojection

• How to map from a cylinder to a planar image?

– Apply camera projection matrix

• w = image width, h = image height

X Y

Z

side view

Convert to image coordinates – Convert to image coordinates

• divide by third coordinate (w)

top-down view

image coords

Cylindrical projection

Y Z X

unit cylinder

Levenberg-Marquardt Method

Alignment

• a rotation of the camera is a translation of the cylinder!

cylinder!

 

 













^

^

 

















 



 

































y y x

x

y y

x

y x

y x y

x

x

y x I y

x J I

y x I y

x J I v

u I

I I

I I I

, 2

, ,

) , ( )

, (

) , ( )

, (



 

 ^{x,y x,y x,y}

LucasKanadeStep

void LucasKanadeStep(CByteImage& img1, CByteImage& img2, float t[2]) { // Transform the image

Translation(img2 img2t t);

Translation(img2, img2t, t);

// Compute the gradients and summed error by comparing img1 and img2t double A[2][2] b[2];

double A[2][2], b[2];

for (int y = 1; y < height-1; y++) { // ignore borders for (int x = 1; x < width-1; x++) {

// If both have full alphas, then compute and accumulate the error // If both have full alphas, then compute and accumulate the error double e = img2t.Pixel(x, y, k) - img1.Pixel (x, y, k);

// Accumulate the matrix entries

double gx = 0.5*(img2t.Pixel(x+1, y, k) - img2t.Pixel(x-1, y, k));g ( g ( , y, ) g ( , y, ));

double gy = 0.5*(img2t.Pixel(x, y+1, k) - img2t.Pixel(x, y-1, k));

A[0][0] += gx*gx; A[0][1] += gx*gy;[ ][ ] g g ; [ ][ ] g gy;

A[1][0] += gx*gy; A[1][1] += gy*gy;

b[0] += e*gx; b[1] += e*gy;

} }

LucasKanadeStep (cont.)

// Solve for the update At=b and update the vector double det 1 0 / (A[0][0]*A[1][1] A[1][0]*A[1][0]);

double det = 1.0 / (A[0][0]*A[1][1] - A[1][0]*A[1][0]);

t[0] += (A[1][1]*b[0] - A[1][0]*b[1]) * det;

t[1] += (A[0][0]*b[1] - A[1][0]*b[0]) * det;

t[1] += (A[0][0] b[1] - A[1][0] b[0]) det;

}

PyramidLucasKanade

void PyramidalLucasKanade(CByteImage& img1, CByteImage& img2, float t[2], int nLevels, int nLucasKanadeSteps)

{{

CBytePyramid p1(img1); // Form the two pyramids CBytePyramid p2(img2);

// Process in a coarse-to-fine hierarchy for (int l = nLevels-1; l >= 0; l--)

{{

t[0] /= (1 << l); // scale the t vector t[1] /= (1 << l);

CByteImage& i1 = p1[l];y g p [ ];

CByteImage& i2 = p2[l];

for (int k = 0; k < nLucasKanadeSteps; k++)( ; p ; ) LucasKanadeStep(i1, i2, t);

t[0] *= (1 << l); // restore the full scaling t[1] *= (1 << l);

} }

Gaussian pyramid

2D Motion models

• translation: x’ = x + t x = (x,y)

i ’ R t

• rotation: x’ = R x + t

• similarity: x’ = s R x + t

• affine: x’ = A x + t

• perspective: x’  H x x = (x y 1)

• perspective: x  H x x = (x,y,1) (x is a homogeneous coordinate)

• These all form a nested group (closed under

• These all form a nested group (closed under composition w/ inv.)

Video matting

alpha matte

Recognising Panoramas

• 1D Rotations ()

• 1D Rotations ()

– Ordering  matching images

Recognising Panoramas

• 1D Rotations ()

• 1D Rotations ()

– Ordering  matching images

Recognising Panoramas

• 1D Rotations ()

• 1D Rotations ()

– Ordering  matching images

Recognising Panoramas

• 1D Rotations ()

• 1D Rotations ()

– Ordering  matching images

• 2D Rotations (q, f)

– Ordering  matching imagesOrdering  matching images

Recognising Panoramas

• 1D Rotations ()

• 1D Rotations ()

– Ordering  matching images

• 2D Rotations (q, f)

– Ordering  matching imagesOrdering  matching images

Recognising Panoramas

• 1D Rotations ()

• 1D Rotations ()

– Ordering  matching images

• 2D Rotations (q, f)

– Ordering  matching imagesOrdering  matching images

Probabilistic model for verification

• Compare probability that this set of RANSAC inliers/outliers was generated by a

inliers/outliers was generated by a correct/false image match

• Choosing values for p₁, p₀ and p_min

Recognising Panoramas

Overview

• SIFT Feature Matching

I M hi

• Image Matching

• Bundle Adjustment

• Multi-band Blending

Nearest Neighbour Matching

• Find k-NN for each feature

k b f l i i ( k 4)

– k  number of overlapping images (we use k = 4)

• Use k-d tree

– k-d tree recursively bi-partitions data at mean in the dimension of maximum variance

A hb f d O l

– Approximate nearest neighbours found in O(nlogn)

Overview

• SIFT Feature Matching

I M hi

• Image Matching

– For each image, use RANSAC to select inlier features

f 6 i ith t f t t h

from 6 images with most feature matches

• Bundle Adjustment

• Multi-band Blending

Finding the panoramas

Finding the panoramas

Finding the panoramas

Finding the panoramas

Overview

• SIFT Feature Matching

I M hi

• Image Matching

• Bundle Adjustment

• Multi-band Blending

Homography for Rotation

• Parameterise each camera by rotation and focal length

focal length

• This gives pairwise homographies

Error function

• Sum of squared projection errors

– n = #imagesn #images

– I(i) = set of image matches to image i

– F(i, j) = set of feature matches between images i,jF(i, j) set of feature matches between images i,j – r_ij^k = residual of k^th feature match between images

i,j,j

• Robust error function

• Robust error function

Overview

• SIFT Feature Matching

I M hi

• Image Matching

• Bundle Adjustment

• Multi-band Blending

Multi-band Blending

• Burt & Adelson 1983

Bl d f b d 

– Blend frequency bands over range  

2-band Blending

Low frequency ( > 2 pixels)

High frequency ( < 2 pixels)

Linear Blending

2-band Blending

Results

Distortion

No distortion Pin cushion Barrel

• Radial distortion of the image

Ca sed b imperfect lenses – Caused by imperfect lenses

– Deviations are most noticeable for rays that pass through the edge of the lens

through the edge of the lens

Radial correction

• Correct for “bending” in wide field of view llenses