%pylab inline
rcParams['figure.figsize'] = (10, 4) #wide graphs by default
from __future__ import print_function
from __future__ import division

from scipy.signal import convolve2d

mat = zeros((100,100))
mat[50,50] = 1
imshow(mat, cmap=cm.gray, interpolation='nearest')

kernel = ones((20,20))/20

from scipy.signal import freqz
ir = ones(20)/20
w, h = freqz(ir)
plot(abs(h))

c = convolve2d(mat, kernel)
imshow(c, cmap= cm.gray, interpolation='nearest')
colorbar()

mat = zeros((100,100))
mat[50,50] = 1
kernel = zeros((20, 20))
kernel[10,  :] = 1
kernel[:, 10] = 1
imshow(kernel, cmap= cm.gray, interpolation='nearest')

c = convolve2d(mat, kernel)
imshow(c, cmap= cm.gray, interpolation='nearest')
colorbar()

a = random.rand(64,64)
snow = where (a > 0.997, 1, 0)
imshow(snow, cmap= cm.gray, interpolation='nearest')

c = convolve2d(snow, kernel)
imshow(c, cmap= cm.gray, interpolation='nearest')
colorbar()

from scipy.ndimage.filters import convolve as convolve2
kernel_ghost = zeros((31, 31))
kernel_ghost[15,15] = 1
kernel_ghost[30,30] = 0.5
imshow(kernel_ghost, cmap = cm.gray, interpolation='nearest')

from scipy.misc import lena
i = lena()
imshow(i, cmap=cm.gray)

i.dtype

ghost = convolve2(i[:,:], kernel_ghost)
imshow(ghost, cmap = cm.gray, interpolation='nearest')
colorbar()

kernel_stamp = zeros((30,30))
kernel_stamp[arange(30), arange(30)] = 1.0
imshow(kernel_stamp, cmap=cm.gray, interpolation='nearest')

subplot(131)
src_image = i[:,:]
imshow(src_image, cmap=cm.gray, interpolation='nearest')

subplot(132)
kernel_stamp = zeros((30,30))
kernel_stamp[arange(30), arange(30)] = 1.0
stamp = convolve2(src_image, kernel_stamp)
imshow(stamp, cmap=cm.gray, interpolation='nearest')

subplot(133)
kernel_stamp = zeros((30,30))
kernel_stamp[15,:] = 1.0
kernel_stamp[:,15] = 1.0
stamp = convolve2(src_image, kernel_stamp)
imshow(stamp, cmap=cm.gray, interpolation='nearest')
gcf().set_figwidth(16)

import scipy.ndimage as ndimage

single_dot = zeros((101, 101))
single_dot[50, 50] = 1.0
gauss_kernel = ndimage.gaussian_filter(single_dot, 5)
imshow(gauss_kernel)
colorbar()

from mpl_toolkits.mplot3d import Axes3D

fig = figure()
ax = Axes3D(fig)
x, y = mgrid[0:101, 0:101]
ax.plot_surface(x,y,gauss_kernel)
fig.add_axes(ax)

plot(abs(fft.rfft(gauss_kernel[:, 50])), 'o-')

plot(abs(fft.rfft(gauss_kernel[:, 50], n=512)))

twinx()
plot(angle(fft.rfft(gauss_kernel[:, 50], n=512)), 'r')

gauss_blur = convolve2(src_image, gauss_kernel)
imshow(gauss_blur, cmap=cm.gray, interpolation= 'nearest')

single_dot = zeros((101, 101))
single_dot[50, 50] = 1.0
gauss_kernel = ndimage.gaussian_filter(single_dot, 20)
imshow(gauss_kernel, interpolation='nearest')
colorbar()

gauss_blur = convolve2(src_image, gauss_kernel)
imshow(gauss_blur, cmap=cm.gray, interpolation= 'nearest')

plot(abs(fft.rfft(gauss_kernel[:, 50], n= 512)))

import scipy.stats as stats
gauss = stats.norm.pdf(arange(17), loc=8, scale=5)
plot(gauss)


kernel = zeros((17,17))
kernel[8,:] = gauss
imshow(kernel, cmap = cm.gray, interpolation='nearest')

plot(abs(rfft(gauss, 512)))

imshow(kernel, interpolation='nearest', cmap=cm.gray)
figure()
imshow(gauss_kernel, interpolation='nearest', cmap=cm.gray)

c = convolve2(src_image, kernel)
subplot(131)
imshow(c, cmap=cm.gray, interpolation='nearest')

c = convolve2(src_image, gauss_kernel)
subplot(132)
imshow(c, cmap=cm.gray, interpolation='nearest')

subplot(133)
imshow(ndimage.gaussian_filter(src_image, 5), cmap=cm.gray, interpolation='nearest')

gcf().set_figwidth(16)

pyt = imread('maxresdefault.jpg')
imshow(pyt)
gcf().set_figheight(8)

from scipy.misc import imresize

resized = imresize(pyt, 0.08, interp='nearest')

imshow(resized, interpolation='nearest')

resized_soft = imresize(ndimage.gaussian_filter(pyt, 4.0), 0.08, interp='nearest')

imshow(resized_soft, interpolation='nearest')

imshow(ndimage.gaussian_filter(resized, 0.4), interpolation='nearest')

single_dot = zeros((101, 101))
single_dot[50, 50] = 1.0
gauss_kernel = ndimage.uniform_filter(single_dot, 20)
imshow(gauss_kernel)
colorbar()

src_image_eyes = src_image[150:300, 100:400]
subplot(121)
imshow(ndimage.gaussian_filter(src_image_eyes, 5), cmap=cm.gray, interpolation='nearest')
subplot(122)
imshow(ndimage.uniform_filter(src_image_eyes, 13), cmap= cm.gray, interpolation='nearest');
gcf().set_figwidth(16)

def imshow2(im):
    imshow(im, cmap=cm.gray, interpolation='nearest')

c3 = ndimage.median_filter(src_image_eyes, size=15)
imshow2(c3)

noise = where(random.random(src_image_eyes.shape) > 0.9, 0.5, 0)
imshow2(noise)

src_image_eyes.max(), src_image_eyes.dtype

noisy_image = (src_image_eyes/255.0) +noise
imshow2(noisy_image)
colorbar()

subplot(131)
imshow2(ndimage.gaussian_filter(noisy_image, 5))
#colorbar()
subplot(132)
imshow2(ndimage.gaussian_filter(noisy_image, 3))
subplot(133)
imshow2(ndimage.gaussian_filter(noisy_image, 1))
gcf().set_figwidth(16)

subplot(311)
imshow2(ndimage.median_filter(noisy_image, 3))
colorbar()

subplot(312)
imshow2(ndimage.median_filter(noisy_image, 10))
colorbar()

subplot(313)
imshow2(ndimage.median_filter(noisy_image, 4))
colorbar()
gcf().set_figheight(8)

title('Sharpen')
subplot(131)
imshow2(src_image)

subplot(132)
c = ndimage.gaussian_filter(src_image, 1)
c = src_image -c*0.5
imshow2(c)

subplot(133)
c = ndimage.gaussian_filter(src_image, 1)
c = src_image -c
imshow2(c)

gauss = stats.norm.pdf(arange(33), loc=16, scale=2.5)
plot(gauss)

sharpen = zeros(33)
sharpen[16] = 1
sharpen -= gauss
plot(sharpen);

plot(abs(rfft(sharpen)));

log_kernel = ndimage.gaussian_laplace(mat, 5)
imshow(log_kernel, interpolation='nearest')
colorbar()

plot(log_kernel[:,50])

plot(abs(fft.rfft(log_kernel[:,50])))

def imshow2(image):
    imshow(image, interpolation='nearest', cmap=cm.gray)

src_image = zeros((100,100))
src_image[25:-25, 25:-25] = 1.0
noise =  where(random.random(src_image.shape) > 0.98, 0.5, 0)
src_image += noise
imshow2(src_image)
colorbar()

structure = zeros((20,20))
structure[range(20), range(20)] = 1.0 #
#structure[range(20), 19 - arange(20)] = 1.0 #
imshow2(structure)

eroded = ndimage.binary_erosion(src_image, structure)
subplot(121)
imshow2(src_image)
subplot(122)
imshow2(eroded)
colorbar()

src_image = where(i[:,:] > 255*0.6, 1, 0)
eroded = ndimage.binary_erosion(src_image, structure)
subplot(121)
imshow2(src_image)
subplot(122)
imshow2(eroded)

eroded = ndimage.binary_erosion(eroded, structure)
subplot(121)
imshow2(eroded)
subplot(122)
eroded = ndimage.binary_erosion(eroded, structure)
imshow2(eroded)

structure = zeros((20,20))
structure[range(20), 19 - arange(20)] = 1.0

src_image = where(i[:,:] > 255*0.6, 1, 0)
eroded = ndimage.binary_erosion(src_image, structure)
subplot(121)
imshow2(src_image)
subplot(122)
imshow2(eroded)

src_image = where(i > 255 * 0.6, 1, 0)
dilated = ndimage.binary_dilation(src_image, structure)
subplot(121)
imshow2(src_image)
subplot(122)
imshow2(dilated)

src_image = where(i > 255 * 0.7, 1, 0)
structure = fliplr(structure)
#structure[range(40), 1-arange(40)]
dilated = ndimage.binary_dilation(src_image, structure)
subplot(121)
imshow2(src_image)
subplot(122)
imshow2(dilated)

mgrid[-3:3, 0:6]

src_image = where(i > 255 * 0.7, 1, 0)
g = mgrid[-15:15, -15:15]
structure = where(sum(g**2, axis=0) < 150, 1, 0)
imshow2(structure)

#structure[range(40), 1-arange(40)]
dilated = ndimage.binary_dilation(src_image, structure)
subplot(121)
imshow2(src_image)
subplot(122)
imshow2(dilated)

src_image = zeros((100,100))
src_image[25:-25, 25:-25] = 1
src_image[10:30, 10:30] = 1
imshow2(src_image)

subplot(121)
imshow2(ndimage.binary_opening(src_image, structure))
title("Opening")
subplot(122)
imshow2(ndimage.binary_closing(src_image, structure))
title("Closing")

imshow2(ndimage.binary_erosion(src_image, structure))

g = mgrid[-5:5, -5:5]
structure = where(sum(g**2, axis=0) < 17, 1, 0)
imshow2(structure)

subplot(121)
imshow2(ndimage.binary_opening(src_image, structure))
title("Opening")
subplot(122)
imshow2(ndimage.binary_closing(src_image, structure))
title("Closing")

imshow2(ndimage.binary_dilation(src_image, structure))

src_image = where(i > 255 * 0.7, 1, 0)
subplot(131)
imshow2(src_image)
subplot(132)
imshow2(ndimage.binary_opening(src_image, structure))
subplot(133)
imshow2(ndimage.binary_closing(src_image, structure))
gcf().set_figwidth(10)

ibwlena = lena()
ibwlena.shape

min_filter = ndimage.minimum_filter(ibwlena, (10,10))
imshow(min_filter, cmap=cm.gray)

subplot(121)
min_filter = ndimage.minimum_filter(ibwlena, (64,64))
imshow(min_filter, cmap=cm.gray)

subplot(122)
max_filter = ndimage.maximum_filter(ibwlena, (64,64))
imshow(max_filter, cmap=cm.gray)

size = (128, 12)
series = ndimage.minimum_filter(ndimage.maximum_filter(ibwlena, size), size)
imshow(series, cmap=cm.gray)

narino = imread("sources/antonio+nariño.jpg")
imshow(narino, cmap=cm.gray)

min_filter = ndimage.minimum_filter(narino, (20,20))
max_filter = ndimage.maximum_filter(narino, (20,20))
imshow(max_filter - min_filter, cmap=cm.gray)

min_filter = ndimage.minimum_filter(narino, (64, 64))
max_filter = ndimage.maximum_filter(narino, (64,64))
imshow(max_filter - min_filter, cmap=cm.gray)

import scipy.ndimage.morphology as morphology

mg = morphology.morphological_gradient(narino, (10,10))
imshow(mg, cmap=cm.gray)

mg = morphology.morphological_gradient(narino, (50,50))
imshow(mg, cmap=cm.gray)

wtp = morphology.white_tophat(narino, (4, 4))
imshow(wtp, cmap=cm.gray);

btp = morphology.black_tophat(narino, (4, 4))
imshow(btp, cmap=cm.gray);

imshow(btp + wtp, cmap=cm.gray);

n = 100
l = 256
im = np.zeros((256,256))
points = l*np.random.random((2, n))
im[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1
im = ndimage.gaussian_filter(im, 6)
im = where(im > im.mean(), 1, 0)
imshow2(im)

label_im, nb_labels = ndimage.label(im)
print(nb_labels) # number of regions

imshow(label_im)
colorbar()

blob = where(label_im == 13, 13, 0)
blob += where(label_im == 22, 22, 0)
imshow(blob, vmax= nb_labels)

blobs = where(label_im > 13, 1, 0)
imshow(blobs)

im = np.zeros((256,256))
im[64:-64, 64:-64] = 1
im = ndimage.rotate(im, 15, mode='constant')
im = ndimage.gaussian_filter(im, 8)
imshow2(im)

sx = ndimage.sobel(im, axis=0, mode='constant')
sy = ndimage.sobel(im, axis=1, mode='constant')
subplot(121)
imshow(sx)
colorbar()
subplot(122)
imshow(sy)
gcf().set_figwidth(16)

sob = hypot(sx, sy)
imshow(sob)
colorbar()

sob = sx + sy
imshow(sob)
colorbar()

imshow(sob > 0.33)

src_mat= zeros((3, 3))
src_mat[1,1] = 1
print(ndimage.sobel(src_mat, axis=0))
print(ndimage.sobel(src_mat, axis=1))

sx = ndimage.prewitt(im, axis=0, mode='constant')
sy = ndimage.prewitt(im, axis=1, mode='constant')
subplot(121)
imshow(sx)
colorbar()
subplot(122)
imshow(sy)
gcf().set_figwidth(16)

pre = hypot(sx, sy)
imshow(pre)
colorbar()

src_mat= zeros((5, 5))
src_mat[2,2] = 1
print(ndimage.sobel(src_mat, axis=0))
print(ndimage.sobel(src_mat, axis=1))
print("-----")
src_mat= zeros((5, 5))
src_mat[2,2] = 1
print(ndimage.prewitt(src_mat, axis=0))
print(ndimage.prewitt(src_mat, axis=1))

src_image = i
sx = ndimage.sobel(src_image, axis=0, mode='constant')
sy = ndimage.sobel(src_image, axis=1, mode='constant')
subplot(121)
imshow2(sx)
colorbar()
subplot(122)
imshow2(sy)

sob = hypot(sx, sy)
imshow2(sob)
colorbar()

imshow2(sob[100:300, 300:450])

src_image = i
sx = ndimage.prewitt(src_image, axis=0, mode='constant')
sy = ndimage.prewitt(src_image, axis=1, mode='constant')
subplot(121)
imshow2(sx)
colorbar()
subplot(122)
imshow2(sy)

pre = hypot(sx, sy)
imshow2(pre)
colorbar()

subplot(121)
imshow2(sob[100:300, 300:450])
colorbar()
subplot(122)
imshow2(pre[100:300, 300:450])
gcf().set_figwidth(16)

src_image = sob
filtered = ndimage.median_filter((src_image > 255 * 0.5), 4)
imshow2(filtered)

src_image = sob[100:300, 300:450]
filtered = ndimage.median_filter((src_image > 255 * 0.5), 4)
imshow2(filtered)

def warp(a, b):
    if a < b:
        out = a + b
    else:
        out = a - b
    return out

vfunc = vectorize(warp)

src_image = i
warped = vfunc(src_image, 127)
imshow(warped, cmap=cm.gray)

hist(src_image.flat, bins=255);

hist(warped.flat, bins=255);

# Assumes uint8 data!
def expand(pixel, center, factor):
    centered = float(pixel) - center
    out = (centered*factor) + center
    if out > 255:
        out = 255
    if out < 0:
        out = 0
    return uint8(out)
    
vfunc = vectorize(expand)

expanded = vfunc(i, 128, 2.0)
imshow(expanded, cmap=cm.gray)

hist(expanded.flat, bins=255);

hist(expanded.flat, bins=255);
ylim((0, 3000))

hist(i.flat, bins=255);
xlim((0, 300));

compressed = vfunc(i, 100, 0.02)
imshow(compressed, cmap=cm.gray)
colorbar()

hist(compressed.flat, bins=255);

compressed = vfunc(i, 10, 0.02)
imshow(compressed, cmap=cm.gray)
colorbar()

hist(compressed.flat, bins=255);

from skimage import color
fruits = imread('vitamy fruits.png')
img_hsv = color.rgb2hsv(fruits)
imshow(fruits)

imshow(img_hsv)

imshow(img_hsv[:,:, 0], cmap= cm.hsv)
colorbar()

imshow(img_hsv[:,:, 1], cmap= cm.gray)

imshow(img_hsv[:,:, 2], cmap= cm.gray)

img_hsv[:,:,0].max()

img_hsv.shape

subplot(231)
title("Hue")
imshow(img_hsv[:,:,0])
subplot(232)
title("Saturation")
imshow(img_hsv[:,:,1])
subplot(233)
title("Value")
imshow(img_hsv[:,:,2])

subplot(234)
imshow(fruits)

low_s = img_hsv[:,:,1] * 0.3
processed = img_hsv.copy()
processed[:,:, 1] = low_s
imshow(color.hsv2rgb(processed))

low_s = img_hsv[:,:,2][:] * 0.3 + 0.7
processed = img_hsv.copy()
processed[:,:, 2] = low_s
imshow(color.hsv2rgb(processed))

low_s = img_hsv[:,:,0][:] + 0.3
processed = img_hsv.copy()
processed[:,:, 0] = low_s
imshow(color.hsv2rgb(processed))

color.

FFT = fft.fft2(narino)

FFT_centered = fft.fftshift(FFT)
    
mag_spec = abs(FFT_centered)
phs_spec = angle(FFT_centered)


imshow(mag_spec, cmap=cm.gray)
colorbar()

imshow(log(mag_spec), cmap=cm.gray)
colorbar()

imshow(phs_spec, cmap=cm.gray)

FFT = fft.fft2(ibwlena)

FFT_centered = fft.fftshift(FFT)
    
mag_spec = abs(FFT_centered)
phs_spec = angle(FFT_centered)

imshow(log(mag_spec), cmap=cm.gray)

imshow(phs_spec, cmap=cm.gray)

FFT1 = fft.fft2(sum(fruits[:400,:332,:], axis=2)/3)
FFT2 = fft.fft2(narino[:400,:332])  

mag_spec = abs(FFT1)
phs_spec = angle(FFT2)

new_FFT = mag_spec * exp(1j*phs_spec)
phs_only = fft.ifft2(new_FFT)
imshow(real(phs_only), cmap=cm.gray)

FFT1 = fft.fft2(sum(fruits[:400,:332,:], axis=2)/3)
FFT2 = fft.fft2(narino[:400,:332])   

mag_spec = abs(FFT2)
phs_spec = angle(FFT1)

new_FFT = mag_spec * exp(1j*phs_spec)
phs_only = fft.ifft2(new_FFT)
imshow(real(phs_only), cmap=cm.gray)