Report

CS 691 Computational Photography Instructor: Gianfranco Doretto Frequency Domain Overview • Frequency domain analysis • Sampling and reconstruction Linear image transformations • In analyzing images, it’s often useful to make a change of basis. F Uf transformed image Vectorized image Fourier transform, or Wavelet transform, or Steerable pyramid transform Self-inverting transforms Same basis functions are used for the inverse transform 1 f U F U F U transpose and complex conjugate Salvador Dali “Gala Contemplating the Mediterranean Sea, which at 30 meters becomes the portrait of Abraham Lincoln”, 1976 A nice set of basis Teases away fast vs. slow changes in the image. This change of basis has a special name… Jean Baptiste Joseph Fourier (1768-1830) manner in which the author arrives at these had crazy idea (1807): ...the equations is not exempt of difficulties and...his Any univariate function can analysis to integrate them still leaves something to be be rewritten as a weighted desired on the score of generality and even rigour. sum of sines and cosines of different frequencies. • Don’t believe it? – Neither did Lagrange, Laplace, Poisson and other big wigs – Not translated into English until 1878! • But it’s (mostly) true! – called Fourier Series – there are some Laplace Lagrange Legendre A sum of sines Our building block: Asin(x Add enough of them to get any signal f(x) you want! Frequency Spectra • example : g(t) = sin(2πf t) + (1/3)sin(2π(3f) t) = + Frequency Spectra Frequency Spectra = = + Frequency Spectra = = + Frequency Spectra = = + Frequency Spectra = = + Frequency Spectra = = + Frequency Spectra 1 = A sin(2 kt ) k 1 k Fourier basis for image analysis Intensity Image Frequency Spectra http://sharp.bu.edu/~slehar/fourier/fourier.html#filtering Signals can be composed + = http://sharp.bu.edu/~slehar/fourier/fourier.html#filtering More: http://www.cs.unm.edu/~brayer/vision/fourier.html Fourier Transform • Fourier transform stores the magnitude and phase at each frequency – Magnitude encodes how much signal there is at a particular frequency – Phase encodes spatial information (indirectly) – For mathematical convenience, this is often notated in terms of real and complex numbers Amplitude: A R( ) I ( ) Euler’s formula: 2 2 I ( ) Phase: tan R( ) 1 Computing the Fourier Transform Continuous Discrete k = -N/2..N/2 Fast Fourier Transform (FFT): NlogN Fourier Transform pairs Phase and Magnitude • Fourier transform of a real function is complex – difficult to plot, visualize – instead, we can think of the phase and magnitude of the transform • Phase is the phase of the complex transform • Magnitude is the magnitude of the complex transform • Curious fact – all natural images have about the same magnitude transform – hence, phase seems to matter, but magnitude largely doesn’t • Demonstration – Take two pictures, swap the phase transforms, compute the inverse - what does the result look like? This is the magnitude transform of the cheetah pic This is the phase transform of the cheetah pic This is the magnitude transform of the zebra pic This is the phase transform of the zebra pic Reconstructio n with zebra phase, cheetah magnitude Reconstruction with cheetah phase, zebra magnitude The Convolution Theorem • The Fourier transform of the convolution of two functions is the product of their Fourier transforms • The inverse Fourier transform of the product of two Fourier transforms is the convolution of the two inverse Fourier transforms • Convolution in spatial domain is equivalent to multiplication in frequency domain! Properties of Fourier Transforms • Linearity • Fourier transform of a real signal is symmetric about the origin • The energy of the signal is the same as the energy of its Fourier transform See Szeliski Book (3.4) Filtering in spatial domain 1 0 -1 2 0 -2 1 0 -1 * = Filtering in frequency domain FFT FFT = Inverse FFT Play with FFT in Matlab • Filtering with fft im = ... % “im” should be a gray-scale floating point image [imh, imw] = size(im); fftsize = 1024; % should be order of 2 (for speed) and include padding im_fft = fft2(im, fftsize, fftsize); % 1) fft im with padding hs = 50; % filter half-size fil = fspecial('gaussian', hs*2+1, 10); fil_fft = fft2(fil, fftsize, fftsize); % 2) fft fil, pad to same size as image im_fil_fft = im_fft .* fil_fft; % 3) multiply fft images im_fil = ifft2(im_fil_fft); % 4) inverse fft2 im_fil = im_fil(1+hs:size(im,1)+hs, 1+hs:size(im, 2)+hs); % 5) remove padding • Displaying with fft figure(1), imagesc(log(abs(fftshift(im_fft)))), axis image, colormap jet Convolution versus FFT • 1-d FFT: O(NlogN) computation time, where N is number of samples. • 2-d FFT: 2N(NlogN), where N is number of pixels on a side • Convolution: K N2, where K is number of samples in kernel • Say N=210, K=100. 2-d FFT: 20 220, while convolution gives 100 220 Why is the Fourier domain particularly useful? • It tells us the effect of linear convolutions. • There is a fast algorithm for performing the DFT, allowing for efficient signal filtering. • The Fourier domain offers an alternative domain for understanding and manipulating the image. coefficient Analysis of our simple filters 1.0 0 Pixel offset original M 1 F [m] f [k ]e Filtered (no change) km i M k 0 1 1.0 0 constant coefficient Analysis of our simple filters 1.0 0 Pixel offset original shifted M 1 F [m] f [k ]e k 0 e i m M km i M Constant magnitude, 1.0 linearly shifted phase 0 coefficient Analysis of our simple filters 0.3 0 Pixel offset original blurred M 1 F [m] f [k ]e km i M k 0 1 m 1 2 cos 3 M 1.0 0 Low-pass filter Analysis of our simple filters 2.0 0.33 0 0 sharpened original M 1 F [m] f [k ]e km i M high-pass filter 2.3 k 0 1 m 2 1 2 cos 3 M 1.0 0 Playing with the DFT of an image Can change spectrum, then reconstruct Low and High Pass filtering Why does the Gaussian give a nice smooth image, but the square filter give edgy artifacts? Gaussian Box filter Gaussian Filter Box Filter Overview • Frequency domain analysis • Sampling and reconstruction Sampling and Reconstruction Sampled representations • How to store and compute with continuous functions? • Common scheme for representation: samples [FvDFH fig.14.14b / Wolberg] – write down the function’s values at many points Reconstruction • Making samples back into a continuous function [FvDFH fig.14.14b / Wolberg] – for output (need realizable method) – for analysis or processing (need mathematical method) – amounts to “guessing” what the function did in between 1D Example: Audio low high frequencies Sampling in digital audio • Recording: sound to analog to samples to disc • Playback: disc to samples to analog to sound again – how can we be sure we are filling in the gaps correctly? Sampling and Reconstruction • Simple example: a sign wave Undersampling • What if we “missed” things between the samples? • Simple example: undersampling a sine wave – unsurprising result: information is lost Undersampling • What if we “missed” things between the samples? • Simple example: undersampling a sine wave – unsurprising result: information is lost – surprising result: indistinguishable from lower frequency Undersampling • What if we “missed” things between the samples? • Simple example: undersampling a sine wave – unsurprising result: information is lost – surprising result: indistinguishable from lower frequency – also was always indistinguishable from higher frequencies – aliasing: signals “traveling in disguise” as other frequencies Aliasing in video Aliasing in images What’s happening? Input signal: Plot as image: x = 0:.05:5; imagesc(sin((2.^x).*x)) Alias! Not enough samples Fourier Transform of a Sampled Signal Sampling period T Sampling frequency F = 1/T T 2T 3T … -F 0 F … Fourier Transform of a Sampled Signal Sampling period T Sampling frequency F = 1/T T 2T 3T … -2F -F 0 F 2F … Nyquist-Shannon Sampling Theorem • When sampling a signal at discrete intervals, the sampling frequency must be F > 2 x fmax • fmax = max frequency of the input signal • This will allow to reconstruct the original perfectly from the sampled version • Do not know fmax or cannot sample at that rate? Prefiltering! We lose something but still better than aliasing!!!!! v v v good bad Image halfsizing This image is too big to fit on the screen. How can we reduce it? How to generate a halfsized version? Image sub-sampling 1/8 1/4 Throw away every other row and column to create a 1/2 size image - called image sub-sampling Image sub-sampling 1/2 1/4 Aliasing! What do we do? (2x zoom) 1/8 (4x zoom) Gaussian (lowpass) pre-filtering G 1/8 G 1/4 Gaussian 1/2 Solution: filter the image, then subsample • Filter size should double for each ½ size reduction. Why? Subsampling with Gaussian pre-filtering Gaussian 1/2 G 1/4 G 1/8 Compare with... 1/2 1/4 (2x zoom) 1/8 (4x zoom) Algorithm for downsampling by factor of 2 1. Start with image(h, w) 2. Apply low-pass filter im_blur = imfilter(image, fspecial(‘gaussian’, 7, 1)) 3. Sample every other pixel im_small = im_blur(1:2:end, 1:2:end); Slide Credits • This set of sides also contains contributions kindly made available by the following authors – Alexei Efros – Frédo Durand – Bill Freeman – Steve Seitz – Derek Hoiem – Steve Marschner