1. Basic examples using the tools in pyrtools

1.1. Content

  • Load an image, and downsample to a size appropriate for the machine speed

  • Synthetic images (‘ramp’, ‘impulse’, etc.)

  • Point operations (lookup tables)

  • histogram Modification/matching

  • Convolution routines

  • Compare speed of convolution/downsampling routines

  • Handling of left and top boundaries (‘reflect1’, ‘reflect2’, ‘repeat’, ‘extend’, ‘zero’, ‘circular’, ‘dont-compute’)

for multi-scale pyramids see TUTORIALS/02_pyramids.ipynb
for multi-scale steerable pyramids see TUTORIALS/03_steerable_pyramids.ipynb
[1]:

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

import pyrtools as pt

%load_ext autoreload
%autoreload 2

1.2. Load an image, and downsample to a size appropriate for the machine speed

[2]:

oim = plt.imread('../DATA/Einstein.pgm').astype(float)

[3]:

imSubSample = 1
im = pt.blurDn(oim, n_levels=imSubSample, filt='qmf9')

[4]:

# imshow: 3 types of automatic graylevel scaling,
# 2 types of automatic sizing, with or without title and Range information.
# to see the documentation for a function, run:`? pt.imshow`

pt.imshow([oim,im], title=['original', 'subsampled'], vrange='auto2', col_wrap=2);
../_images/tutorials_01_tools_5_0.png

1.3. Statistics:

[5]:

pt.image_stats(im);

Image statistics:
  Range: [10.126404, 245.754672]
  Mean: 116.771655,  Stdev: 37.977892,  Kurtosis: 3.085949

1.4. Synthetic images

[6]:

# pick some parameters
size      = 256
direction = 2 * np.pi * np.random.rand(1)
slope     = 10 * np.random.rand(1) - 5
intercept = 10 * np.random.rand(1) - 5
origin    = np.round(size * np.random.rand(2,1)).astype(int)
exponent  = 0.8 + np.random.rand(1)
amplitude = 1 + 5 * np.random.rand(1)
phase     = 2 * np.pi * np.random.rand(1)
period    = 20
twidth    = 7

[7]:

pt.imshow(pt.synthetic_images.ramp(size, direction, slope, intercept, origin));
../_images/tutorials_01_tools_10_0.png 
[8]:

pt.imshow(pt.synthetic_images.impulse(size, origin, amplitude));
../_images/tutorials_01_tools_11_0.png 
[9]:

pt.imshow(pt.synthetic_images.polar_radius(size, exponent, origin));
../_images/tutorials_01_tools_12_0.png 
[10]:

pt.imshow(pt.synthetic_images.polar_angle(size, direction));
../_images/tutorials_01_tools_13_0.png 
[11]:

pt.imshow(pt.synthetic_images.disk(size, size/4, origin, twidth));
../_images/tutorials_01_tools_14_0.png 
[12]:

pt.imshow(pt.synthetic_images.gaussian(size, (size/6)**2, origin, 'norm'));
../_images/tutorials_01_tools_15_0.png 
[13]:

# TODO fix normalization - range
g = pt.synthetic_images.gaussian(size, (size/6)**2, origin, 'norm')
g.min(), g.max()

[13]:

(9.402979666533684e-14, 8.742642137616321e-05)

[14]:

pt.imshow(pt.synthetic_images.gaussian(size, ((size/8)**2,(size/3)**2), origin, amplitude));
../_images/tutorials_01_tools_17_0.png 
[15]:

# cov = (size * np.random.uniform(-1,1,(2,2)))
# cov = cov.dot(cov.T)
# print(np.round(cov))
cov = np.array((( 12571., -25233.),
                (-25233.,  52488.)))

pt.imshow(pt.synthetic_images.gaussian(size, cov, origin, amplitude));
../_images/tutorials_01_tools_18_0.png 
[16]:

fig=pt.imshow(pt.synthetic_images.zone_plate(size, amplitude, phase));
fig.savefig('zoneplate.png')
../_images/tutorials_01_tools_19_0.png 
[17]:

pt.imshow(pt.synthetic_images.angular_sine(size, 3, amplitude, phase, origin));
../_images/tutorials_01_tools_20_0.png 
[18]:

pt.imshow(pt.synthetic_images.sine(size, period, direction, amplitude=amplitude, phase=phase, origin=origin));
../_images/tutorials_01_tools_21_0.png 
[19]:

pt.imshow(pt.synthetic_images.sine(20, frequency=[1,2]));
../_images/tutorials_01_tools_22_0.png 
[20]:

pt.imshow(pt.synthetic_images.square_wave(size, period, direction, amplitude, phase=phase, origin=origin, twidth=twidth));
../_images/tutorials_01_tools_23_0.png 
[21]:

pt.imshow(pt.synthetic_images.pink_noise(size, exponent));
../_images/tutorials_01_tools_24_0.png 
[22]:

# Plotting complex valued images

sp = np.fft.fftshift(np.fft.fft2(im, norm='ortho'))

pt.imshow([im, sp], plot_complex='logpolar');
../_images/tutorials_01_tools_25_0.png

1.5. Point operations (lookup tables):

[23]:

Xtbl,Ytbl = pt.rcosFn(width=20, position=25, values=(-1, 1))
plt.plot(Xtbl,Ytbl)
plt.show()
../_images/tutorials_01_tools_27_0.png 
[24]:

pt.imshow(pt.pointOp(pt.synthetic_images.polar_radius(size,1,[70,30]), Ytbl, Xtbl[0], Xtbl[1]-Xtbl[0], 0));
../_images/tutorials_01_tools_28_0.png

1.6. Convolution routines: Compare speed of convolution/downsampling routines

[25]:

k = 5
size = 2 ** 9
noise = np.random.rand(size,size)
filt  = np.random.rand(k,k)

[26]:

%%time
res1 = pt.corrDn(noise, np.flipud(np.fliplr(filt)), 'reflect1', step=[2, 2])

CPU times: user 4.54 ms, sys: 8.26 ms, total: 12.8 ms
Wall time: 4.15 ms

[27]:

%%time
ires = pt.rconv2(noise,filt) # this is using scipy.signal.convolve
res2 = ires[0:size:2,0:size:2]

CPU times: user 48.9 ms, sys: 18.1 ms, total: 67 ms
Wall time: 22.3 ms

[28]:

res1 = pt.corrDn(noise, np.flipud(np.fliplr(filt)), 'reflect1', step=[2, 2])
ires = pt.rconv2(noise,filt) # this is using scipy.signal.convolve
res2 = ires[0:size:2,0:size:2]
pt.image_compare(res1, res2);

Difference statistics:
  Range: [0, 0]
  Mean: -0.000000,  Stdev (rmse): 0.000000,  SNR (dB): 294.068252

[29]:

res3 = pt.corrDn(noise, np.flipud(np.fliplr(filt)), 'reflect1', step=[1, 1])
pt.image_compare(res3, ires);

Difference statistics:
  Range: [0, 0]
  Mean: -0.000000,  Stdev (rmse): 0.000000,  SNR (dB): 294.006866

[30]:

im = plt.imread('../DATA/Einstein.pgm').astype(float)
# im = pt.synthetic_images.impulse(256, origin=None)

[31]:

binom5 = pt.binomial_filter(5)
# construct a separable 2D filter
lo_filt = 2*binom5*binom5.T

[32]:

blurred = pt.rconv2(im, lo_filt)
pt.imshow([im, blurred], title=['im', 'blurred']);
../_images/tutorials_01_tools_37_0.png 
[33]:

# much faster implementation:
b = pt.corrDn(im, np.flipud(np.fliplr(lo_filt)), 'reflect1', step=[1, 1])
pt.image_compare(blurred, b);

# each dimension separately
bx  = pt.corrDn(im, binom5.T, 'reflect1')
bxy = pt.corrDn(bx, binom5, 'reflect1')
bxy *= 2
pt.image_compare(bxy, b);

Difference statistics:
  Range: [0, 0]
  Mean: 0.000000,  Stdev (rmse): 0.000000,  SNR (dB): 305.282881
Difference statistics:
  Range: [0, 0]
  Mean: 0.000000,  Stdev (rmse): 0.000000,  SNR (dB): inf

[34]:

# blur and downsample in a single step
blurred1 = pt.corrDn(image = im, filt = lo_filt, step = (2,2))

bx  = pt.corrDn(im, binom5.T, 'reflect1', step=[1, 2])
bxy = pt.corrDn(bx, binom5, 'reflect1', step=[2, 1])
bxy *= 2

pt.imshow([blurred1, bxy], zoom=2);
pt.image_compare(blurred1, bxy);

Difference statistics:
  Range: [0, 0]
  Mean: 0.000000,  Stdev (rmse): 0.000000,  SNR (dB): inf
../_images/tutorials_01_tools_39_1.png 
[35]:

# NOTE: here is helper to decide on subsampling depending on machine speed

import time
t = time.time()
pt.corrDn(oim, filt=np.ones((2,2)) / 4, edge_type='reflect1', step=(2, 2), start=(0, 0), stop=None)
elapsed = time.time() - t
imSubSample = min(max(np.floor(np.log2(elapsed)/2+3),0),2)
imSubSample

[35]:

0

1.7. Handling of left and top boundaries:

[36]:

fsz = [9, 9]
fmid = np.ceil((fsz[0]+1)/2)
imsz = (16, 16)

# pick one:
im = np.eye(imsz[0])
im = pt.synthetic_images.ramp(imsz, np.pi/6)
im = pt.synthetic_images.square_wave(imsz, 6, np.pi/6)

# pick one:
edges='reflect1'
edges='reflect2'
edges='repeat'
edges='extend'
edges='zero'
edges='circular'
edges='dont-compute'

filt = pt.synthetic_images.impulse(fsz,[0, 0])
pt.imshow(pt.corrDn(im,filt,edges), zoom=10, title='Edges: ' + str(edges))
# TODO
plt.plot([0,0,imsz[1],imsz[1],0]+fmid-1.5,
         [0,imsz[0],imsz[0],0,0]+fmid-1.5, lw=.5)
plt.show()
../_images/tutorials_01_tools_42_0.png