Azimuth image tutorial

%auto # We need functions from numpy and scipy/pylab in Sage. import numpy as np, pylab as pl # read in the image the Data directory of this worksheet cat = pl.imread(DATA+'leo.png') # here we doing a Fourier transform, which Numpy supports jazz_cat = np.fft.fft2(cat) # inverse Fourier transform to try to get the cat back cf = np.fft.ifft2(jazz_cat) # here we construct a list of graphic objects - for now matrix plots mat_cats = [matrix_plot(im) for im in [cat,real_part(jazz_cat), real_part(cf), cat - real_part(cf)]] ga = graphics_array(mat_cats,2,2) html("<h3>Cat photo and how it travels to the spectral domain and back again.</h3>") html("<h4>Original photo: size %.2e and variance %.2f</h4>" %(cat.size,cat.var())) ga.show( axes=False, figsize=(12,6), aspect_ratio = 1) 

Cat photo and how it travels to the spectral domain and back again.
Original photo: size 2.14e+06 and variance 0.06

Cat photo and how it travels to the spectral domain and back again.
Original photo: size 2.14e+06 and variance 0.06

timeit('np.fft.fft2(cat)') 

5 loops, best of 3: 1.16 s per loop

5 loops, best of 3: 1.16 s per loop

# answers to hw q.2 cat.shape, prod(cat.shape) == cat.size,cat.min() 

((840, 849, 3), True, 0.0)

((840, 849, 3), True, 0.0)

# siam 1.2.4, fig 1.2 but color cat matrix_plot(.1*(2*cat -1) ^4+0.4) 

# one way to read an image from the web directly import urllib f=urllib.urlopen("http://imgs.xkcd.com/comics/literally.png") img = f.read() 

# convert the image to scalar form - one scalar value/pixel, mean along two axes bw_cat = pl.mean(cat,2) # make a copy of the black and white - scalar - image copy_cat = np.copy(bw_cat) noisy_cat = copy_cat + 0.1*np.random.normal(size=copy_cat.shape) # time-series in Sage "understands" numpy n-arrays cs = finance.TimeSeries(bw_cat, False) g = graphics_array([matrix_plot(bw_cat), matrix_plot(noisy_cat), cs.plot_histogram()],2,3) show(g, axes=False, figsize=(8,6), aspect_ratio = 1) html("<h3>A copy cat in black and white and a tad noisy cat</h3>") 

A copy cat in black and white and a tad noisy cat

A copy cat in black and white and a tad noisy cat

 

import scipy.signal 

# from the "scipy cookbook" at http://www.scipy.org/Cookbook def gauss_kern(size, sizey=None): """ Returns a normalized 2D gauss kernel array for convolutions """ size = int(size) if not sizey: sizey = size else: sizey = int(sizey) x, y = np.mgrid[-size:size+1, -sizey:sizey+1] g = np.exp(-(x**2/float(size)+y**2/float(sizey))) return g / g.sum() def blur_image(im, n, ny=None) : """ blurs the image by convolving with a gaussian kernel of typical size n. The optional keyword argument ny allows for a different size in the y direction. """ g = gauss_kern(n, sizey=ny) improc = scipy.signal.convolve(im,g, mode='valid') return(improc) 

X, Y = np.mgrid[-70:70, -70:70] Z = np.cos((X**2+Y**2)/200.) + np.random.normal(size=X.shape) matrix_plot(Z) 

matrix_plot(blur_image(noisy_cat,3)) 

time _ = blur_image(noisy_cat,3) 

Time: CPU 0.85 s, Wall: 0.84 s

Time: CPU 0.85 s, Wall: 0.84 s

a= animate([matrix_plot(blur_image(noisy_cat,w)) for w in srange(1,22,2)],figsize=8) 

a.show(delay=17) 

%hide def sgolay2d ( z, window_size, order, derivative=None): """ """ # number of terms in the polynomial expression n_terms = ( order + 1 ) * ( order + 2) / 2.0 if window_size % 2 == 0: raise ValueError('window_size must be odd') if window_size**2 < n_terms: raise ValueError('order is too high for the window size') half_size = window_size // 2 # exponents of the polynomial. # p(x,y) = a0 + a1*x + a2*y + a3*x^2 + a4*y^2 + a5*x*y + ... # this line gives a list of two item tuple. Each tuple contains # the exponents of the k-th term. First element of tuple is for x # second element for y. # Ex. exps = [(0,0), (1,0), (0,1), (2,0), (1,1), (0,2), ...] exps = [ (k-n, n) for k in range(order+1) for n in range(k+1) ] # coordinates of points ind = np.arange(-half_size, half_size+1, dtype=np.float64) dx = np.repeat( ind, window_size ) dy = np.tile( ind, [window_size, 1]).reshape(window_size**2, ) # build matrix of system of equation A = np.empty( (window_size**2, len(exps)) ) for i, exp in enumerate( exps ): A[:,i] = (dx**exp[0]) * (dy**exp[1]) # pad input array with appropriate values at the four borders new_shape = z.shape[0] + 2*half_size, z.shape[1] + 2*half_size Z = np.zeros( (new_shape) ) # top band band = z[0, :] Z[:half_size, half_size:-half_size] = band - np.abs( np.flipud( z[1:half_size+1, :] ) - band ) # bottom band band = z[-1, :] Z[-half_size:, half_size:-half_size] = band + np.abs( np.flipud( z[-half_size-1:-1, :] ) -band ) # left band band = np.tile( z[:,0].reshape(-1,1), [1,half_size]) Z[half_size:-half_size, :half_size] = band - np.abs( np.fliplr( z[:, 1:half_size+1] ) - band ) # right band band = np.tile( z[:,-1].reshape(-1,1), [1,half_size] ) Z[half_size:-half_size, -half_size:] = band + np.abs( np.fliplr( z[:, -half_size-1:-1] ) - band ) # central band Z[half_size:-half_size, half_size:-half_size] = z # top left corner band = z[0,0] Z[:half_size,:half_size] = band - np.abs( np.flipud(np.fliplr(z[1:half_size+1,1:half_size+1]) ) - band ) # bottom right corner band = z[-1,-1] Z[-half_size:,-half_size:] = band + np.abs( np.flipud(np.fliplr(z[-half_size-1:-1,-half_size-1:-1]) ) - band ) # top right corner band = Z[half_size,-half_size:] Z[:half_size,-half_size:] = band - np.abs( np.flipud(Z[half_size+1:2*half_size+1,-half_size:]) - band ) # bottom left corner band = Z[-half_size:,half_size].reshape(-1,1) Z[-half_size:,:half_size] = band - np.abs( np.fliplr(Z[-half_size:, half_size+1:2*half_size+1]) - band ) # solve system and convolve if derivative == None: m = np.linalg.pinv(A)[0].reshape((window_size, -1)) return scipy.signal.fftconvolve(Z, m, mode='valid') elif derivative == 'col': c = np.linalg.pinv(A)[1].reshape((window_size, -1)) return scipy.signal.fftconvolve(Z, -c, mode='valid') elif derivative == 'row': r = np.linalg.pinv(A)[2].reshape((window_size, -1)) return scipy.signal.fftconvolve(Z, -r, mode='valid') elif derivative == 'both': c = np.linalg.pinv(A)[1].reshape((window_size, -1)) r = np.linalg.pinv(A)[2].reshape((window_size, -1)) return scipy.signal.fftconvolve(Z, -r, mode='valid'), scipy.signal.fftconvolve(Z, -c, mode='valid') 

# create some sample twoD data x = np.linspace(-3,3,100) y = np.linspace(-3,3,100) X, Y = np.meshgrid(x,y) Z = np.exp( -(X**2+Y**2)) # add noise Zn = Z + np.random.normal( 0, 0.2, Z.shape ) # filter it Zf = sgolay2d( Zn, window_size=29, order=4) # plot original, original plus noise, and filtered mlist = [matrix_plot(n) for n in [Z,Zn,Zf]] 

ga = graphics_array(mlist) ga.show() 

@interact def _(w=(1,101,2), order=(1,20,1)): show(matrix_plot(sgolay2d(noisy_cat,window_size = w, order=order)),figsize=8) 

a= animate([matrix_plot(sgolay2d(noisy_cat,window_size = w, order=2)) for w in srange(3,77,2)],figsize=8) 

# 60 fps or 17 ms a.show(delay=17) 

# making an published interact of SG # We need functions from numpy and scipy/pylab in Sage. import numpy as np, pylab as pl import scipy.signal # read in the image the Data directory of this worksheet cat = pl.imread(DATA+'leo.png') # convert the image to scalar form - one scalar value/pixel, mean along two axes bw_cat = pl.mean(cat,2) # make a copy of the black and white - scalar - image copy_cat = np.copy(bw_cat) noisy_cat = copy_cat + 0.1*np.random.normal(size=copy_cat.shape) def sgolay2d ( z, window_size, order, derivative=None): """ """ # number of terms in the polynomial expression n_terms = ( order + 1 ) * ( order + 2) / 2.0 if window_size % 2 == 0: raise ValueError('window_size must be odd') if window_size**2 < n_terms: raise ValueError('order is too high for the window size') half_size = window_size // 2 # exponents of the polynomial. # p(x,y) = a0 + a1*x + a2*y + a3*x^2 + a4*y^2 + a5*x*y + ... # this line gives a list of two item tuple. Each tuple contains # the exponents of the k-th term. First element of tuple is for x # second element for y. # Ex. exps = [(0,0), (1,0), (0,1), (2,0), (1,1), (0,2), ...] exps = [ (k-n, n) for k in range(order+1) for n in range(k+1) ] # coordinates of points ind = np.arange(-half_size, half_size+1, dtype=np.float64) dx = np.repeat( ind, window_size ) dy = np.tile( ind, [window_size, 1]).reshape(window_size**2, ) # build matrix of system of equation A = np.empty( (window_size**2, len(exps)) ) for i, exp in enumerate( exps ): A[:,i] = (dx**exp[0]) * (dy**exp[1]) # pad input array with appropriate values at the four borders new_shape = z.shape[0] + 2*half_size, z.shape[1] + 2*half_size Z = np.zeros( (new_shape) ) # top band band = z[0, :] Z[:half_size, half_size:-half_size] = band - np.abs( np.flipud( z[1:half_size+1, :] ) - band ) # bottom band band = z[-1, :] Z[-half_size:, half_size:-half_size] = band + np.abs( np.flipud( z[-half_size-1:-1, :] ) -band ) # left band band = np.tile( z[:,0].reshape(-1,1), [1,half_size]) Z[half_size:-half_size, :half_size] = band - np.abs( np.fliplr( z[:, 1:half_size+1] ) - band ) # right band band = np.tile( z[:,-1].reshape(-1,1), [1,half_size] ) Z[half_size:-half_size, -half_size:] = band + np.abs( np.fliplr( z[:, -half_size-1:-1] ) - band ) # central band Z[half_size:-half_size, half_size:-half_size] = z # top left corner band = z[0,0] Z[:half_size,:half_size] = band - np.abs( np.flipud(np.fliplr(z[1:half_size+1,1:half_size+1]) ) - band ) # bottom right corner band = z[-1,-1] Z[-half_size:,-half_size:] = band + np.abs( np.flipud(np.fliplr(z[-half_size-1:-1,-half_size-1:-1]) ) - band ) # top right corner band = Z[half_size,-half_size:] Z[:half_size,-half_size:] = band - np.abs( np.flipud(Z[half_size+1:2*half_size+1,-half_size:]) - band ) # bottom left corner band = Z[-half_size:,half_size].reshape(-1,1) Z[-half_size:,:half_size] = band - np.abs( np.fliplr(Z[-half_size:, half_size+1:2*half_size+1]) - band ) # solve system and convolve if derivative == None: m = np.linalg.pinv(A)[0].reshape((window_size, -1)) return scipy.signal.fftconvolve(Z, m, mode='valid') elif derivative == 'col': c = np.linalg.pinv(A)[1].reshape((window_size, -1)) return scipy.signal.fftconvolve(Z, -c, mode='valid') elif derivative == 'row': r = np.linalg.pinv(A)[2].reshape((window_size, -1)) return scipy.signal.fftconvolve(Z, -r, mode='valid') elif derivative == 'both': c = np.linalg.pinv(A)[1].reshape((window_size, -1)) r = np.linalg.pinv(A)[2].reshape((window_size, -1)) return scipy.signal.fftconvolve(Z, -r, mode='valid'), scipy.signal.fftconvolve(Z, -c, mode='valid') @interact def _(w=(1,101,2), order=(1,20,1)): show(matrix_plot(sgolay2d(noisy_cat,window_size = w, order=order)),figsize=8)